      program main

c*********************************************************************72
c
cc MAIN is the main program for TOMS644_PRB.
c
c  Discussion:
c
c    TOMS644_PRB tests the TOMS644 library.
c
c  Modified:
c
c    13 February 2010
c
      implicit none

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS644_PRB:'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test the TOMS644 library.'
c
c  Single precision
c
      call cqcbh ( )
      call cqcbi ( )
      call cqcbj ( )
      call cqcbk ( )
      call cqcby ( )
      call cqcai ( )
c
c  Double precision
c
      call zqcbh ( )
      call zqcbi ( )
      call zqcbj ( )
      call zqcbk ( )
      call zqcby ( )
      call zqcai ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS644_PRB:'
      write ( *, '(a)' ) '  Normal end of execution.'

      stop
      end
      subroutine CQCBH

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C     CQCBH IS A QUICK CHECK ROUTINE FOR THE COMPLEX H BESSEL FUNCTIONS
C     GENERATED BY SUBROUTINE CBESH.
C
C     CQCBH GENERATES SEQUENCES OF H BESSEL FUNCTIONS FOR KIND 2 FROM
C     CBESH AND CHECKS THEM AGAINST ANALYTIC CONTINUATION FORMULAS
C     IN THE (Z,FNU) SPACE:
C
C     KODE = 1 TESTS (ANALYTIC CONTINUATION FORMULAE, I**2 = -1):
C
C     H(FNU,2,Z)=-EXP(I*PI*FNU)*H(FNU,1,-Z),       -PI.LT.ARG(Z).LE.0
C
C               = 2*COS(PI*FNU)*H(FNU,2,-Z) + EXP(I*PI*FNU)*H(FNU,1,-Z),
C
C                                                   0.LT.ARG(Z).LE.PI
C
C     KODE = 2 TESTS FOR KINDS 1 AND 2:
C
C            EXP(-I*Z)*H(FNU,1,Z) = [EXP(-I*Z)*H(FNU,1,Z)]
C
C            EXP( I*Z)*H(FNU,2,Z) = [EXP( I*Z)*H(FNU,2,Z)]
C
C     WHERE THE LEFT SIDE IS COMPUTED WITH KODE = 1 AND THE RIGHT SIDE
C     WITH KODE = 2.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
      COMPLEX CW, CI, U, V, W, Y, Z, ZN, CSGN
      REAL AA, AB, AER, ALIM, ATOL, AV, CT, DIG, ERR, XX, YY,
     * ELIM, EPS, ER, ERTOL, FNU, FNUL, PI, R, RL,
     * RM, R1M4, R1M5, R2, ST, T, TOL, TS, XNU, R1MACH, SLAK, FILM
      INTEGER I, ICASE, IERR, IHP, IL, IR, IRB, IT, ITL, K, KODE, KK,
     *K1, K2, LFLG, LUN, MFLG, M, N, NU, NZ1, NZ2, NZ3, I1MACH, KEPS,
     *MQC, NL, NUL, KDO
      DIMENSION T(20), AER(20), XNU(20), U(20), V(20), W(20), Y(20),
     *KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = AMAX1(R1M4,1.0E-18)
      AA = -ALOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN0(IABS(K1),IABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + AMAX1(-AB,-41.45E0)
      DIG = AMIN1(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0
      SLAK = AMAX1(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = AMIN1(RM,200.0E0)
      RM = AMAX1(RM,RL+10.0E0)
      R2 = AMIN1(FNUL,RM)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE H BESSEL FUNCTIONS FROM CBES
     *H'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6E12.4/)
      ATOL = 100.0E0*TOL
      PI = 4.0E0*ATAN(1.0E0)
      CI = CMPLX(0.0E0,1.0E0)
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.1E0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.1E0
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM=FLOAT(IL-1)
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*FLOAT(-IL+2*K-1)/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 170 KODE=1,2
        DO 160 N=1,NL
          DO 150 NU=1,NUL
            FNU = XNU(NU)
            DO 140 ICASE=1,3
              IRB = MIN0(ICASE,2)
              DO 130 IR=IRB,3
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (EPS*FLOAT(3-IR)+2.0E0*FLOAT(IR-1))/2.0E0
                GO TO 80
   60           CONTINUE
                R = (2.0E0*FLOAT(3-IR)+R2*FLOAT(IR-1))/2.0E0
                GO TO 80
   70           CONTINUE
                IF (R2.GE.RM) GO TO 140
                R = (R2*FLOAT(3-IR)+RM*FLOAT(IR-1))/2.0E0
   80           CONTINUE
                DO 120 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0E0
                  IF (ABS(ST).LT.ATOL) ST = 0.0E0
                  XX = R*CT
                  YY = R*ST
                  Z = CMPLX(XX,YY)
                  IF (KODE.EQ.1) THEN
                    M=2
                    CALL CBESH(Z,FNU,KODE,M,N,Y,NZ1,IERR)
                    IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120
                    IF (ST.LT.0.0E0 .OR. (ST.EQ.0.0E0.AND.CT.GT.0.0E0))
     *              THEN
                      IHP = 1
                      ZN=-Z
                      M=1
                      CALL CBESH(ZN,FNU,KODE,M,N,W,NZ2,IERR)
                      IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    ELSE
                      IHP = 2
                      ZN=-Z
                      M=2
                      CALL CBESH(ZN,FNU,KODE,M,N,W,NZ3,IERR)
                      IF (IERR.NE.0.OR.NZ3.NE.0) GO TO 120
                      M=1
                      CALL CBESH(ZN,FNU,KODE,M,N,V,NZ2,IERR)
                      IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    ENDIF
                    AB=AMOD(FNU,2.0E0)*PI
                    CSGN = CMPLX(COS(AB),SIN(AB))
                    MFLG = 0
                    DO 100 I=1,N
                      AB = FNU+FLOAT(I-1)
                      AA = MAX(0.5E0,AB)
                      IF (IHP.EQ.1) THEN
                        V(I) = -CSGN*W(I)
                        CW = Y(I) - V(I)
                      ELSE
                        V(I) = 2.0E0*REAL(CSGN)*W(I) + CSGN*V(I)
                        CW = Y(I) - V(I)
                      ENDIF
                      AV = CABS(Y(I))
                      ER = CABS(CW)
                      IF (YY.EQ.0.0E0) THEN
                        IF (ABS(XX).LT.AA) ER = ER/AV
                      ELSE
                        ER = ER/AV
                      ENDIF
                      AER(I) = ER
                      IF (ER.GT.ERTOL) MFLG = 1
                      CSGN = -CSGN
  100               CONTINUE
                  ELSE
                    M=1
                    KK=1
                    CALL CBESH(Z,FNU,KK,M,N,U,NZ1,IERR)
                    IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120
                    CALL CBESH(Z,FNU,KODE,M,N,V,NZ2,IERR)
                    IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    M=2
                    KK=1
                    CALL CBESH(Z,FNU,KK,M,N,W,NZ1,IERR)
                    IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120
                    CALL CBESH(Z,FNU,KODE,M,N,Y,NZ2,IERR)
                    IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    ZN=CI*Z
                    ZN=EXP(ZN)
                    MFLG = 0
                    DO 105 I=1,N
                      AB = FNU+FLOAT(I-1)
                      AA = MAX(0.5E0,AB)
                      CW = U(I)/ZN-V(I)
                      ER = CABS(CW)
                      AV = CABS(V(I))
                      IF (YY.EQ.0.0E0) THEN
                        IF (ABS(XX).LT.AA) ER = ER/AV
                      ELSE
                        ER = ER/AV
                      ENDIF
                      ERR = ER
                      IF (ER.GT.ERTOL) MFLG = 1
                      CW=ZN*W(I)-Y(I)
                      ER = CABS(CW)
                      AV = CABS(Y(I))
                      IF (YY.EQ.0.0E0) THEN
                        IF (ABS(XX).LT.AA) ER = ER/AV
                      ELSE
                        ER = ER/AV
                      ENDIF
                      IF (ER.GT.ERTOL) MFLG = 1
                      AER(I) = ER+ERR
  105               CONTINUE
                  ENDIF
                  IF (MFLG.EQ.0) GO TO 120
                  IF (LFLG.EQ.1) GO TO 110
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,V(1),Y(1)')
                  LFLG = 1
  110             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, ICASE
99992             FORMAT (5I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) Z, FNU, V(1), Y(1)
99991             FORMAT (7E12.4)
  120           CONTINUE
  130         CONTINUE
  140       CONTINUE
  150     CONTINUE
  160   CONTINUE
  170 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine CQCBI

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C     CQCBI IS A QUICK CHECK ROUTINE FOR THE COMPLEX I BESSEL FUNCTION
C     GENERATED BY SUBROUTINE CBESI.
C
C     CQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM
C     CBESI AND CBESK AND CHECKS THE WRONSKIAN EVALUATION
C
C           I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z
C
C     IN THE RIGHT HALF PLANE AND A MODIFIED FORM
C
C          I(FNU+1,Z)*K(FNU,ZR) - I(FNU,Z)*K(FNU+1,ZR) = C/Z
C
C     IN THE LEFT HALF PLANE WHERE ZR=-Z AND C=EXP(I*FNU*SGN) WITH
C     SGN=+1 FOR IM(Z).GE.0 AND SGN=-1 FOR IM(Z).LT.0.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
      COMPLEX  CONE, CSGN, CC, CV, CW, CY, W, Y, Z, ZR
      REAL AA, AB, AER, ALIM, ARG, ATOL, DIG, ELIM, EPS,
     * ER, ERTOL, FFNU, FNU, FNUL, GNU, HPI, PI, R, RL, R1M4, R1M5,
     * R2, RM, T, TOL, XX, YY, R1MACH, SLAK, TS, ST, CT, FILM, XNU
      INTEGER I, ICASE, IFNU, IL, IPRNT, IR, IT, ITL, IRB, K, KK, KODE,
     * K1, K2, LFLG, LUN, MFLG, N, NU, NZ, N1, NUL, IERR,
     * MQC, NL, KEPS, KDO, I1MACH
      DIMENSION T(20), AER(20), Y(20), W(20), XNU(20), KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = MAX(R1M4,1.0E-18)
      AA = -ALOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN(ABS(K1),ABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + MAX(-AB,-41.45E0)
      DIG = MIN(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0
      SLAK = MAX(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = MIN(RM,200.0E0)
      RM = MAX(RM,RL+10.0E0)
      R2 = MIN(FNUL,RM)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM CBESI
     *'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6E12.4/)
      CONE = CMPLX(1.0E0,0.0E0)
      ATOL = 100.0E0*TOL
      HPI = 2.0E0*ATAN(1.0E0)
      PI = HPI + HPI
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.1E0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.1E0
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM=FLOAT(IL-1)
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*FLOAT(-IL+2*K-1)/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 200 KODE=1,2
        DO 190 N=1,NL
          N1 = N + 1
          DO 180 NU=1,NUL
            FNU = XNU(NU)
            IFNU = INT(FNU)
            FFNU = FNU - FLOAT(IFNU)
            ARG = PI*FFNU
            CSGN = CMPLX(COS(ARG),SIN(ARG))
            IF (MOD(IFNU,2).EQ.1) CSGN = -CSGN
            DO 170 ICASE=1,3
              IRB = MIN(2,ICASE)
              DO 160 IR=IRB,4
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (0.2E0*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0
                GO TO 80
   60           CONTINUE
                R = (2.0E0*FLOAT(4-IR)+R2*FLOAT(IR-1))/3.0E0
                GO TO 80
   70           CONTINUE
                IF (R2.GE.RM) GO TO 170
                R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0
   80           CONTINUE
                DO 150 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0E0
                  IF (ABS(ST).LT.ATOL) ST = 0.0E0
                  XX = R*CT
                  YY = R*ST
                  Z = CMPLX(XX,YY)
                  IF (CT.GE.0.0E0) THEN
C-----------------------------------------------------------------------
C     WRONSKIAN CHECKS IN THE RIGHT HALF PLANE
C-----------------------------------------------------------------------
                    CALL CBESI(Z, FNU, KODE, N1, Y, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CALL CBESK(Z, FNU, KODE, N1, W, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS
C     ON KODE=2
C-----------------------------------------------------------------------
                    CV = CONE/Z
                    IF (KODE.EQ.2) THEN
                      CW = CMPLX(COS(YY),SIN(YY))
                      CV = CW*CV
                    ENDIF
                    CC = CONE
                  ELSE
C-----------------------------------------------------------------------
C     WRONSKIAN CHECKS IN THE LEFT HALF PLANE
C-----------------------------------------------------------------------
                    ZR = -Z
                    CALL CBESI(Z, FNU, KODE, N1, Y, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CALL CBESK(ZR, FNU, KODE, N1, W, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CV = CSGN
                    IF (YY.LT.0.0E0) THEN
                      CV = CONJG(CV)
                    ENDIF
                    CV = CV/Z
                    IF (KODE.EQ.2) THEN
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS
C     ON KODE=2. SCALE FACTOR = EXP(-I*YY) FOR RE(Z).LT.0
C-----------------------------------------------------------------------
                      CW = CMPLX(COS(YY),-SIN(YY))
                      CV = CV*CW
                    ENDIF
                    CC = -CONE
                  ENDIF
                  MFLG = 0
                  KK=0
                  DO 130 I=1,N
                    CW = W(I)*Y(I+1)
                    CY = CC*W(I+1)*Y(I)
                    CY = CY + CW - CV
                    ER = CABS(CY)/CABS(CV)
                    AER(I) = ER
                    IF (ER.GT.ERTOL) THEN
                      IF(KK.EQ.0) THEN
                        MFLG = 1
                        KK=I
                      ENDIF
                    ENDIF
                    IF (CT.LT.0.0E0) CV = -CV
  130             CONTINUE
                  IF (MFLG.EQ.0) GO TO 150
                  IF (LFLG.EQ.1) GO TO 140
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK)        KK=INDEX O
     *F FIRST NON-ZERO PAIR'/)
                  LFLG = 1
  140             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, ICASE, KK
99992             FORMAT (6I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) Z, FNU, Y(KK)
99991             FORMAT (6E12.4)
  150           CONTINUE
  160         CONTINUE
  170       CONTINUE
  180     CONTINUE
  190   CONTINUE
  200 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      IF (MQC.EQ.1) then
        return
      end if
C-----------------------------------------------------------------------
C     CHECKS NEAR UNDERFLOW LIMITS ON SERIES(I=1) AND UNIFORM
C     ASYMPTOTIC EXPANSION(I=2)
C-----------------------------------------------------------------------
      write ( *,99989)
99989 FORMAT (/' CHECKS NEAR UNDERFLOW AND OVERFLOW LIMITS'/)
      Z = CMPLX(1.4E0,1.4E0)
      IPRNT = 0
      DO 280 I=1,2
        FNU = 10.2E0
        KODE = 1
        N = 20
  230   CONTINUE
        CALL CBESI(Z, FNU, KODE, N, Y, NZ, IERR)
        IF (NZ.NE.0) GO TO 240
        FNU = FNU + 5.0E0
        GO TO 230
  240   CONTINUE
        IF (NZ.LT.10) GO TO 250
        FNU = FNU - 1.0E0
        GO TO 230
  250   CONTINUE
        CALL CBESK(Z, FNU, KODE, 2, W, NZ, IERR)
        CV = CONE/Z
        CY = W(1)*Y(2)
        CW = W(2)*Y(1)
        CW = CW + CY - CV
        ER = CABS(CW)/CABS(CV)
        IF (ER.LT.ERTOL) GO TO 270
        IF (IPRNT.EQ.1) GO TO 260
        write ( *,99988)
99988   FORMAT (/' OUTPUT FORMAT'/' ERROR,Z,FNU,KODE,N'/)
        IPRNT = 1
  260   CONTINUE
        write ( *,99987) ER, Z, FNU, KODE, N
99987   FORMAT (4E12.4, 2I5)
  270   CONTINUE
        XX = RL + RL
        Z = CMPLX(XX,0.0E0)
  280 CONTINUE
C-----------------------------------------------------------------------
C     CHECK NEAR OVERFLOW LIMITS
C-----------------------------------------------------------------------
      Z = CMPLX(ELIM,0.0E0)
      FNU = 0.0
  290 CONTINUE
      CALL CBESK(Z, FNU, KODE, N, Y, NZ, IERR)
      IF (NZ.LT.10) GO TO 300
      IF (NZ.EQ.N) FNU = FNU + 3.0E0
      FNU = FNU + 2.0E0
      GO TO 290
  300 CONTINUE
      GNU = FNU + FLOAT(N-2)
      CALL CBESI(Z, GNU, KODE, 2, W, NZ, IERR)
      CV = CONE/Z
      CY = Y(N-1)*W(2)
      CW = Y(N)*W(1)
      CW = CW + CY - CV
      ER = CABS(CW)/CABS(CV)
      IF (ER.LT.ERTOL) GO TO 320
      IF (IPRNT.EQ.1) GO TO 310
      write ( *,99988)
      IPRNT = 1
  310 CONTINUE
      write ( *,99987) ER, Z, FNU, KODE, N
  320 CONTINUE
      IF (IPRNT.EQ.0) write ( *,99990)
      return
      END
      subroutine CQCBJ

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C     CQCBJ IS A QUICK CHECK ROUTINE FOR THE COMPLEX J BESSEL FUNCTION
C     GENERATED BY SUBROUTINE CBESJ.
C
C     CQCBJ GENERATES SEQUENCES OF J AND H BESSEL FUNCTIONS FROM CBESJ
C     AND CBESH AND CHECKS THE WRONSKIANS
C
C     J(FNU,Z)*H(FNU+1,1,Z)-J(FNU+1,Z)*H(FNU,1,Z)=2/(PI*I*Z)   Y.GE.0
C
C     J(FNU,Z)*H(FNU+1,2,Z)-J(FNU+1,Z)*H(FNU,2,Z)=-2/(PI*I*Z)  Y.LT.0
C
C     IN THEIR RESPECTIVE HALF PLANES.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
      COMPLEX Z, WR, CJ, CH, CON, T1, T2, CER
      REAL AA, AB, AER, ALIM, DIG, ELIM, EPS, ER, ERTOL, FNU, FNUL,
     * GNU, HPI, R, RL, RM, R1M4, R1M5, R2, T, TOL, XNU, XX, YY,
     * R1MACH, SLAK, FILM, ST, TS, CT, SGN
      INTEGER I, ICASE, IL, IR, IRB, IT, ITL, K, KK, KODE, K1, K2,
     * LFLG, LUN, M, MFLG, N, NU, NZJ, NZH, IERRJ, IERRH, I1MACH, KEPS,
     * KDO, NL, NUL, MQC
      DIMENSION T(20), AER(20), XNU(20), CJ(20), CH(20), KEPS(20),
     * KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = MAX(R1M4,1.0E-18)
      AA = -LOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN(ABS(K1),ABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + MAX(-AB,-41.45E0)
      DIG = MIN(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-LOG10(TOL)-7.0E0)/11.0E0
      SLAK = MAX(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = MIN(RM,200.0E0)
      RM = MAX(RM,RL+10.0E0)
      R2 = MIN(FNUL,RM)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM CBESJ
     *'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6E12.4/)
      ATOL=100.0E0*TOL
      HPI = 2.0E0*ATAN(1.0E0)
      PI = HPI + HPI
      CON=CMPLX(0.0,-1.0/HPI)
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.1E0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.1E0
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM=FLOAT(IL-1)
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*FLOAT(-IL+2*K-1)/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 260 KODE=1,2
        DO 250 N=1,NL
          NP=N+1
          DO 240 NU=1,NUL
            FNU = XNU(NU)
            GNU = FNU + FLOAT(N-1) + 1.0E0
            GNU=SQRT(GNU)
            GNU=MIN(GNU,0.5*RL)
            DO 230 ICASE=1,3
              IRB = MIN(2,ICASE)
              DO 220 IR=IRB,4
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (GNU*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0
                GO TO 80
   60           CONTINUE
                R = (2.0E0*FLOAT(4-IR)+R2*FLOAT(IR-1))/3.0E0
                GO TO 80
   70           CONTINUE
                IF (R2.GE.RM) GO TO 230
                R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0
   80           CONTINUE
                DO 210 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0E0
                  IF (ABS(ST).LT.ATOL) ST = 0.0E0
                  Z = CMPLX(R*CT,R*ST)
                  XX = REAL(Z)
                  YY = AIMAG(Z)
                  IF(XX.EQ.0.0.AND.YY.EQ.0.0) GO TO 210
                  WR=CON/Z
                  M=1
                  IF(YY.LT.0.0) THEN
                    M=2
                    WR=-WR
                  ENDIF
                  CALL CBESJ(Z,FNU,KODE,NP,CJ,NZJ,IERRJ)
                  CALL CBESH(Z,FNU,KODE,M,NP,CH,NZH,IERRH)
                  IF(NZJ.NE.0.OR.NZH.NE.0) GO TO 210
                  IF(IERRJ.NE.0.OR.IERRH.NE.0) GO TO 210
                  IF(KODE.EQ.2) THEN
                    SGN=3.0-2.0*FLOAT(M)
                    WR=WR*CMPLX(COS(XX),-SGN*SIN(XX))
                  ENDIF
                  KK = 0
                  MFLG = 0
                  DO 190 I=1,N
                    T1=CJ(I)*CH(I+1)
                    T2=CJ(I+1)*CH(I)
                    CER=T1-T2-WR
                    ER=CABS(CER)/CABS(WR)
                    IF (ER.GT.ERTOL) THEN
                      IF(MFLG.EQ.0) THEN
                        MFLG = 1
                        KK=I
                      ENDIF
                    ENDIF
                    AER(I)=ER
  190             CONTINUE
                  IF (MFLG.EQ.0) GO TO 210
                  IF (LFLG.EQ.1) GO TO 200
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZJ,NZH,ICASE
     *')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,CJ(KK),CH(KK), KK=INDEX
     * OF FIRST NON-ZERO CJ,CH PAIR'/)
                  LFLG = 1
  200             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, NZJ, NZH, ICASE
99992             FORMAT (8I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) Z, FNU, CJ(KK), CH(KK)
99991             FORMAT (9E12.4)
  210           CONTINUE
  220         CONTINUE
  230       CONTINUE
  240     CONTINUE
  250   CONTINUE
  260 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine CQCBK

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C     CQCBK IS A QUICK CHECK ROUTINE FOR THE COMPLEX K BESSEL FUNCTION
C     GENERATED BY SUBROUTINE CBESK.
C
C     CQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM
C     CBESI AND CBESK AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION
C
C           I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z
C
C     IN THE RIGHT HALF PLANE AND THE ANALYTIC CONTINUATION FORMULA
C     FOR H(FNU,2,Z) IN TERMS OF THE K FUNCTION
C
C           K(FNU,Z) = C3*H(FNU,2,ZR) + C4*H(FNU,1,ZR)    IM(Z).GE.0
C
C                    = CONJG(K(FNU,CONJG(Z)))             IM(Z).LT.0
C
C     IN THE LEFT HALF PLANE WHERE C3=C1*CONJG(C2)*C5, C4 = C2*C5
C     C1=2*COS(PI*FNU),   C2=EXP(PI*FNU*I/2),   C5 =-PI*I/2   AND
C     ZR = Z*EXP(-3*PI*I/2) = Z*I
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
      COMPLEX CONE, CSGN, CV, CW, CY, C1, C2, C3, C4, V, W, Y, Z, ZR,
     * ZZ, CIP, COE
      REAL AA, AB, AER, ALIM, ARG, ATOL, DIG, ELIM, EPS, ER,
     * ERTOL, FFNU, FNU, FNUL, HPI, PI, R, RL, RM, R1M4, R1M5, R2,
     * T, TOL, XNU, XX, R1MACH, FILM, ST, CT, TS, SLAK, YY
      INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, I4, K, KK, KODE,
     * K1, K2, LFLG, LUN, M, MFLG, N, NU, NZ, N1, IERR, I1MACH,
     * KEPS, KDO, MQC, NL, NUL
      DIMENSION T(20), AER(20), XNU(20), V(20), Y(20), W(20), KEPS(20),
     * KDO(20), CIP(4)

      DATA CIP(1),CIP(2),CIP(3),CIP(4) / (1.0E0,0.0E0),(0.0E0,1.0E0),
     * (-1.0E0,0.0E0),(0.0E0,-1.0E0) /
      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = MAX(R1M4,1.0E-18)
      AA = -ALOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN(ABS(K1),ABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + MAX(-AB,-41.45E0)
      DIG = MIN(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0
      SLAK = MAX(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = MIN(RM,200.0E0)
      RM = MAX(RM,RL+10.0E0)
      R2 = MIN(FNUL,RM)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE K BESSEL FUNCTION FROM CBESK
     *'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6E12.4/)
      CONE = CMPLX(1.0E0,0.0E0)
      ATOL = 100.0E0*TOL
      HPI = 2.0E0*ATAN(1.0E0)
      PI = HPI + HPI
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.1E0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.1E0
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM=FLOAT(IL-1)
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*FLOAT(-IL+2*K-1)/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 200 KODE=1,2
        DO 190 N=1,NL
          N1 = N + 1
          DO 180 NU=1,NUL
            FNU = XNU(NU)
            IFNU = INT(FNU)
            FFNU = FNU - FLOAT(IFNU)
            ARG = HPI*FFNU
            CSGN = CMPLX(COS(ARG),SIN(ARG))
            I4 = MOD(IFNU,4)+1
            CSGN = CSGN*CIP(I4)
            DO 170 ICASE=1,3
              IRB = MIN(2,ICASE)
              DO 160 IR=IRB,4
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (0.2E0*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0
                GO TO 80
   60           CONTINUE
                R = (2.0E0*FLOAT(4-IR)+R2*FLOAT(IR-1))/3.0E0
                GO TO 80
   70           CONTINUE
                IF (R2.GE.RM) GO TO 170
                R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0
   80           CONTINUE
                DO 150 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0E0
                  IF (ABS(ST).LT.ATOL) ST = 0.0E0
                  XX = R*CT
                  YY = R*ST
                  Z = CMPLX(XX,YY)
                  IF (XX.GE.0.0E0) THEN
C-----------------------------------------------------------------------
C     WRONSKIAN CHECKS IN THE RIGHT HALF PLANE
C-----------------------------------------------------------------------
                    CALL CBESI(Z, FNU, KODE, N1, W, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CALL CBESK(Z, FNU, KODE, N1, Y, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS
C     ON KODE=2
C-----------------------------------------------------------------------
                    CV = CONE/Z
                    IF (KODE.EQ.2) THEN
                      CV = CV*CMPLX(COS(YY),SIN(YY))
                    ENDIF
                    MFLG = 0
                    KK=0
                    DO 130 I=1,N
                      CW = W(I)*Y(I+1)
                      CY = W(I+1)*Y(I)
                      CY = CY + CW - CV
                      ER = CABS(CY)/CABS(CV)
                      AER(I) = ER
                      IF (ER.GT.ERTOL) THEN
                        IF(KK.EQ.0) THEN
                          MFLG = 1
                          KK=I
                        ENDIF
                      ENDIF
  130               CONTINUE
                  ELSE
C-----------------------------------------------------------------------
C     ANALYTIC CONTINUATION FORMULA CHECKS FOR LEFT HALF PLANE IN TERMS
C     OF H(FNU,1,Z) AND H(FNU,2,Z)
C-----------------------------------------------------------------------
                    ZZ = Z
                    IF (YY.LT.0.0E0) THEN
                      ZZ=CONJG(Z)
                    ENDIF
                    ZR = ZZ*CMPLX(0.0E0,1.0E0)
                    M=1
                    CALL CBESH(ZR, FNU, KODE, M, N, W, NZ, IERR)
                    IF (IERR.NE.0) GO TO 150
                    M=2
                    CALL CBESH(ZR, FNU, KODE, M, N, V, NZ, IERR)
                    IF (IERR.NE.0) GO TO 150
                    CALL CBESK(Z, FNU, KODE, N, Y, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    COE=CMPLX(0.0E0,-HPI)
                    MFLG = 0
                    KK = 0
                    AA = 2.0E0*COS(PI*FFNU)
                    IF(MOD(IFNU,2).NE.0) AA = -AA
                    C1 = CMPLX(AA,0.0E0)
                    C2 = CSGN
                    DO 135 I=1,N
                      C3 = C1
                      C4 = C2
                      IF (KODE.EQ.2) THEN
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO COEFICIENTS DUE TO SCALING OF H(FNU,1,Z) AND
C     H(FNU,2,Z) FUNCTIONS ON KODE = 2.
C-----------------------------------------------------------------------
                        AB=CABS(V(I))
                        AA=ALOG(AB)+XX+XX
                        IF (AA.GT.ELIM) GO TO 150
                        IF (AA.LT.-ELIM) THEN
                           C3 = CMPLX(0.0E0,0.0E0)
                        ELSE
                          CW = ZZ+ZZ
                          C3 = C3*CEXP(CW)
                        ENDIF
                      ENDIF
                      CY = (C3*CONJG(C2)*V(I)+C4*W(I))*COE
                      IF (YY.LT.0.0E0) THEN
                        CY=CONJG(CY)
                      ENDIF
                      ER = CABS(CY-Y(I))/CABS(Y(I))
                      AER(I) = ER
                      IF (ER.GT.ERTOL) THEN
                        IF(KK.EQ.0) THEN
                          MFLG = 1
                          KK=I
                        ENDIF
                      ENDIF
                      C2 = C2*CMPLX(0.0E0,1.0E0)
                      C1 = -C1
  135               CONTINUE
                  ENDIF
                  IF (MFLG.EQ.0) GO TO 150
                  IF (LFLG.EQ.1) GO TO 140
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK)        KK=INDEX O
     *F FIRST NON-ZERO PAIR'/)
                  LFLG = 1
  140             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, ICASE, KK
99992             FORMAT (6I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) Z, FNU, Y(KK)
99991             FORMAT (6E12.4)
  150           CONTINUE
  160         CONTINUE
  170       CONTINUE
  180     CONTINUE
  190   CONTINUE
  200 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine CQCBY

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C     CQCBY IS A QUICK CHECK ROUTINE FOR THE COMPLEX Y BESSEL FUNCTION
C     GENERATED BY SUBROUTINE CBESY.
C
C     CQCBY GENERATES SEQUENCES OF Y BESSEL FUNCTIONS FROM CBESY AND
C     CBESYH AND COMPARES THEM FOR A VARIETY OF VALUES IN THE (Z,FNU)
C     SPACE. CBESYH IS AN OLD VERSION OF CBESY WHICH COMPUTES THE Y
C     FUNCTION FROM THE H FUNCTIONS OF KINDS 1 AND 2.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
      COMPLEX  CW, CWRK, V, W, Z
      REAL AA, AB, AER, ALIM, ATOL, AV, CABS, DIG, ELIM, EPS, ER,
     * ERTOL, FFNU, FNU, FNUL, PI, R, RL, RM, R1M4, R1M5, R2,
     * T, TOL, XNU, XX, YY, R1MACH, CT, ST, TS, SLAK, FILM
      INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, K, KK,
     * KODE, K1, K2, LFLG, LUN, MFLG, N, NU, NZ1, NZ2, IERR, I1MACH,
     * KEPS, KDO, NL, NUL, MQC
      DIMENSION T(20), AER(20), XNU(20), W(20), V(20),
     * CWRK(20), KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = AMAX1(R1M4,1.0E-18)
      AA = -ALOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN0(IABS(K1),IABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + AMAX1(-AB,-41.45E0)
      DIG = AMIN1(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0
      SLAK = AMAX1(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = AMIN1(RM,200.0E0)
      RM = AMAX1(RM,RL+10.0E0)
      R2 = AMIN1(FNUL,RM)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM CBESY
     * '/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6E12.4/)
      ATOL = 100.0E0*TOL
      PI = 4.0E0*ATAN(1.0E0)
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI/2.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.1E0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.1E0
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM=FLOAT(IL-1)
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*FLOAT(-IL+2*K-1)/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 190 KODE=1,2
        DO 180 N=1,NL
          DO 170 NU=1,NUL
            FNU = XNU(NU)
            IFNU = INT(FNU)
            FFNU = FNU-FLOAT(IFNU)
            DO 160 ICASE=1,3
              IRB = MIN0(2,ICASE)
              DO 150 IR=IRB,4
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (EPS*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0
                GO TO 80
   60           CONTINUE
                R = (2.0E0*FLOAT(4-IR)+R2*FLOAT(IR-1))/3.0E0
                GO TO 80
   70           CONTINUE
                IF (R2.EQ.RM) GO TO 160
                R = (R2*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0
   80           CONTINUE
                DO 140 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0E0
                  IF (ABS(ST).LT.ATOL) ST = 0.0E0
                  Z = CMPLX(R*CT,R*ST)
                  XX = REAL(Z)
                  YY = AIMAG(Z)
                  CALL CBESY(Z, FNU, KODE, N, V, NZ2, CWRK, IERR)
                  IF (NZ2.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL CBESYH(Z, FNU, KODE, N, W, NZ1, CWRK, IERR)
                  IF (NZ1.NE.0.OR.IERR.NE.0) GO TO 140
                  MFLG = 0
                  DO 120 I=1,N
                    AB = FNU+FLOAT(I-1)
                    AA = MAX(0.5E0,AB)
                    CW = W(I) - V(I)
                    AV = CABS(V(I))
                    ER = CABS(CW)
                    IF (AV.NE.0.0E0) THEN
                      IF (YY.EQ.0.0E0) THEN
                        IF (XX.GT.0.0E0) THEN
                          IF (ABS(XX).LT.AA) ER = ER/AV
                        ELSE
                          IF (ABS(FFNU-0.5E0).LT.0.125E0) THEN
                            IF (ABS(XX).LT.AA) ER = ER/AV
                          ELSE
                            ER = ER/AV
                          ENDIF
                        ENDIF
                      ELSE
                        ER = ER/AV
                      ENDIF
                    ENDIF
                    AER(I) = ER
                    IF (ER.GT.ERTOL) MFLG = 1
  120             CONTINUE
                  IF (MFLG.EQ.0) GO TO 140
                  IF (LFLG.EQ.1) GO TO 130
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL = ', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZ1,NZ2,ICASE
     *')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,W(KK),V(KK), KK=INDEX O
     *F FIRST NON-ZERO W,V PAIR'/)
                  LFLG = 1
  130             CONTINUE
                  KK = MAX0(NZ1,NZ2) + 1
                  KK = MIN0(N,KK)
                  write ( *,99992) KODE, N, IR, IT, NZ1, NZ2, ICASE
99992             FORMAT (8I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) Z, FNU, W(KK), V(KK)
99991             FORMAT (7E12.4)
  140           CONTINUE
  150         CONTINUE
  160       CONTINUE
  170     CONTINUE
  180   CONTINUE
  190 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      SUBROUTINE CBESYH(Z, FNU, KODE, N, CY, NZ, CWRK, IERR)

c*********************************************************************72
c
C***BEGIN PROLOGUE  CBESYH
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C***CATEGORY NO.  B5K
C***KEYWORDS  Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
C             BESSEL FUNCTION OF SECOND KIND
C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
C***PURPOSE  TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
C***DESCRIPTION
C
C         ON KODE=1, CBESYH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
C         BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
C         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
C         -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESYH RETURNS THE SCALED
C         FUNCTIONS
C
C         CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
C
C         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
C         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
C         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
C         (REF. 1).
C
C         INPUT
C           Z      - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI
C           FNU    - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0E0
C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
C                    KODE= 1  RETURNS
C                             CY(I)=Y(FNU+I-1,Z), I=1,...,N
C                        = 2  RETURNS
C                             CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
C                             WHERE Y=AIMAG(Z)
C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
C           CWRK   - A COMPLEX WORK VECTOR OF DIMENSION AT LEAST N
C
C         OUTPUT
C           CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
C                    VALUES FOR THE SEQUENCE
C                    CY(I)=Y(FNU+I-1,Z)  OR
C                    CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y))  I=1,...,N
C                    DEPENDING ON KODE.
C           NZ     - NZ=0 , A NORMAL RETURN
C                    NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
C                    UNDERFLOW (GENERALLY ON KODE=2)
C           IERR   - ERROR FLAG
C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
C                    IERR=1, INPUT ERROR   - NO COMPUTATION
C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU+N-1 IS
C                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
C                            ACCURACY
C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
C                            CANCE BY ARGUMENT REDUCTION
C                    IERR=5, ERROR              - NO COMPUTATION,
C                            ALGORITHM TERMINATION CONDITION NOT MET
C
C***LONG DESCRIPTION
C
C         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
C
C              Y(FNU,Z) = 0.5*(H(1,FNU,Z) - H(2,FNU,Z))/I
C
C         WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
C         AND H(2,FNU,Z) ARE CALCULATED IN CBESH.
C
C         FOR NEGATIVE ORDERS,THE FORMULA
C
C              Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
C
C         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
C         INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
C         POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
C         SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
C         NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
C         LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
C         CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
C         WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
C         ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
C
C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
C         IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
C
C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
C         OR -PI/2+P.
C
C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
C                 COMMERCE, 1955.
C
C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
C
C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
C
C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
C                 1018, MAY, 1985
C
C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
C                 MATH. SOFTWARE, 12, NO. 3, SEPTEMBER 1986, PP 265-273.
C
C***ROUTINES CALLED  CBESH,I1MACH,R1MACH
C***END PROLOGUE  CBESYH
C
      COMPLEX CWRK, CY, C1, C2, EX, HCI, Z, ZU, ZV
      REAL ELIM, EY, FNU, R1, R2, TAY, XX, YY, R1MACH, ASCLE, RTOL,
     * ATOL, AA, BB, R1M5, TOL
      INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
      DIMENSION CY(N), CWRK(N)
C***FIRST EXECUTABLE STATEMENT  CBESYH
      XX = REAL(Z)
      YY = AIMAG(Z)
      IERR = 0
      NZ=0
      IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1
      IF (FNU.LT.0.0E0) IERR=1
      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
      IF (N.LT.1) IERR=1
      IF (IERR.NE.0) RETURN
      HCI = CMPLX(0.0E0,0.5E0)
      CALL CBESH(Z, FNU, KODE, 1, N, CY, NZ1, IERR)
      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
      CALL CBESH(Z, FNU, KODE, 2, N, CWRK, NZ2, IERR)
      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
      NZ = MIN0(NZ1,NZ2)
      IF (KODE.EQ.2) GO TO 60
      DO 50 I=1,N
        CY(I) = HCI*(CWRK(I)-CY(I))
   50 CONTINUE
      RETURN
   60 CONTINUE
      TOL = AMAX1(R1MACH(4),1.0E-18)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      K = MIN0(IABS(K1),IABS(K2))
      R1M5 = R1MACH(5)
C-----------------------------------------------------------------------
C     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
C-----------------------------------------------------------------------
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      R1 = COS(XX)
      R2 = SIN(XX)
      EX = CMPLX(R1,R2)
      EY = 0.0E0
      TAY = ABS(YY+YY)
      IF (TAY.LT.ELIM) EY = EXP(-TAY)
      IF (YY.LT.0.0E0) GO TO 90
      C1 = EX*CMPLX(EY,0.0E0)
      C2 = CONJG(EX)
   70 CONTINUE
      NZ = 0
      RTOL = 1.0E0/TOL
      ASCLE = R1MACH(1)*RTOL*1.0E+3
      DO 80 I=1,N
C       CY(I) = HCI*(C2*CWRK(I)-C1*CY(I))
        ZV = CWRK(I)
        AA=REAL(ZV)
        BB=AIMAG(ZV)
        ATOL=1.0E0
        IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 75
          ZV = ZV*CMPLX(RTOL,0.0E0)
          ATOL = TOL
   75   CONTINUE
        ZV = ZV*C2*HCI
        ZV = ZV*CMPLX(ATOL,0.0E0)
        ZU=CY(I)
        AA=REAL(ZU)
        BB=AIMAG(ZU)
        ATOL=1.0E0
        IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 85
          ZU = ZU*CMPLX(RTOL,0.0E0)
          ATOL = TOL
   85   CONTINUE
        ZU = ZU*C1*HCI
        ZU = ZU*CMPLX(ATOL,0.0E0)
        CY(I) = ZV - ZU
        IF (CY(I).EQ.CMPLX(0.0E0,0.0E0) .AND. EY.EQ.0.0E0) NZ = NZ + 1
   80 CONTINUE
      RETURN
   90 CONTINUE
      C1 = EX
      C2 = CONJG(EX)*CMPLX(EY,0.0E0)
      GO TO 70
  170 CONTINUE
      NZ = 0
      RETURN
      END
      subroutine CQCAI

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C     CQCAI IS A QUICK CHECK ROUTINE FOR THE COMPLEX AIRY FUNCTIONS
C     GENERATED BY SUBROUTINES CAIRY AND CBIRY.
C
C     CQCAI GENERATES AIRY FUNCTIONS AND THEIR DERIVATIVES FROM CAIRY
C     AND CBIRY AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION IN THE
C     REGION -PI/3 .LE. ARG(Z) .LE. PI/3:
C
C                 AI(Z)*BI'(Z)-AI'(Z)*BI(Z)=1/PI.
C
C     IN THE REMAINDER OF THE CUT PLANE, THE IDENTITIES
C
C              AI(Z)  = SQRT(-Z)*( J(-1/3,ZR) + J(1/3,ZR) )/3
C
C              AI'(Z) =        Z*( J(-2/3,ZR) - J(2/3,ZR) )/3
C
C       BI(Z)  =   I*SQRT(-Z/3)*( C1*H(1/3,1,ZR) - C2*H(1/3,2,ZR) )/2
C
C       BI'(Z) = I*(-Z)/SQRT(3)*( C2*H(2/3,1,ZR) - C1*H(2/3,2,ZR) )/2
C
C     ARE CHECKED WHERE ZR = (2/3)(-Z)**(3/2) WITH C1 = EXP(PI*I/6),
C     C2 = CONJG(C1) AND I**2 = -1.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
      COMPLEX CA, CAV, CHI, CI, CONA, CONB, CONC, COND, CON1, CON2, CV,
     *  CW, CY, W, Y, YY, Z, ZR, ZW, SC
      REAL AA, AB, ALIM, ATOL, AV, CT, C13, C23, C43, ZX, ZY,
     * DIG, ELIM, EPS, ER, ERTOL, FNUL, FPI, HPI, PI, R, RL, RPI,
     * R1M4, R1M5, SPI, ST, SLAK, T, TOL, RM, FILM, PI3, R1MACH, TS
      INTEGER I, ICASE, ICL, IL, IR, IRSET, IT, ITL, J, JB, JL, K,
     * KEPS, KODE, K1, K2, LFLG, LUN, NZ, IERR, I1MACH, MQC, IRB, KDO
      DIMENSION  ER(5), T(20), Y(20), YY(20), W(20), KDO(20), KEPS(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = R1MACH(4)
      TOL = MAX(R1M4,1.0E-18)
      AA = -ALOG10(R1M4)
      K1 = I1MACH(12)
      K2 = I1MACH(13)
      R1M5 = R1MACH(5)
      K = MIN(ABS(K1),ABS(K2))
      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
      AB = AA*2.303E0
      ALIM = ELIM + MAX(-AB,-41.45E0)
      DIG = MIN(AA,18.0E0)
      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
      RL = 1.2E0*DIG + 3.0E0
      SLAK = 3.0E0+4.0E0*(-ALOG10(TOL)-7.0E0)/11.0E0
      SLAK = MAX(SLAK,3.0E0)
      ERTOL = TOL*10.0E0**SLAK
      RM = 0.5E0*(ALIM + ELIM)
      RM = MIN(RM,200.0E0)
      RM = MAX(RM,RL+10.0E0)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE AIRY FUNCTIONS FROM CAIRY AN
     *D CBIRY'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6E12.4/)
      ATOL = 100.0E0*TOL
      FPI = ATAN(1.0E0)
      HPI = FPI + FPI
      PI = HPI + HPI
      RPI = 1.0E0/PI
      SPI = PI/6.0E0
      CON1 = CMPLX(COS(SPI),SIN(SPI))
      CON2 = CONJG(CON1)
      PI3 = SPI+SPI
      C13 = 1.0E0/3.0E0
      C23 = C13+C13
      C43 = C23+C23
      AV = SQRT(C13)
      CAV = CMPLX(AV,0.0E0)
      CHI = CMPLX(0.0E0,0.5E0)
      CI = CMPLX(0.0E0,1.0E0)
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        ICL=1
        IL=5
        DO 5 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    5   CONTINUE
      ELSE
        ICL=2
        IL=7
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KEPS(2)=1
        KEPS(3)=1
        KEPS(5)=1
        KEPS(6)=1
      ENDIF
      I = 2
      EPS = 0.01E0
      FILM=FLOAT(IL-1)
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF(KDO(K).EQ.0) THEN
          T(I) = PI*FLOAT(-IL+2*K-1)/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 180 ICASE=1,ICL
        DO 170 KODE=1,2
          DO 160 IRSET=1,3
            IRB = MIN(IRSET,2)
            DO 150 IR=IRB,4
              GO TO (40, 50, 60), IRSET
   40         CONTINUE
              R = (0.2E0*FLOAT(4-IR)+2.0E0*FLOAT(IR-1))/3.0E0
              GO TO 70
   50         CONTINUE
              R = (2.0E0*FLOAT(4-IR)+RL*FLOAT(IR-1))/3.0E0
              GO TO 70
   60         CONTINUE
              R = (RL*FLOAT(4-IR)+RM*FLOAT(IR-1))/3.0E0
   70         CONTINUE
              DO 140 IT=1,ITL
                CT = COS(T(IT))
                ST = SIN(T(IT))
                IF (ABS(CT).LT.ATOL) CT = 0.0E0
                IF (ABS(ST).LT.ATOL) ST = 0.0E0
                ZX = R*CT
                ZY = R*ST
                Z = CMPLX(ZX,ZY)
                IF(ABS(T(IT)).LE.PI3) THEN
C-----------------------------------------------------------------------
C     WRONSKIAN CHECK IN -PI/3.LT.ARG(Z).LT.PI/3, TEST #1
C-----------------------------------------------------------------------
                  CALL CAIRY(Z, 0, KODE, Y(1), NZ, IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL CAIRY(Z, 1, KODE, Y(2), NZ, IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL CBIRY(Z, 0, KODE, W(1), IERR)
                  CALL CBIRY(Z, 1, KODE, W(2), IERR)
                  CW = Y(1)*W(2)
                  CY = Y(2)*W(1)
                  CV = CMPLX(RPI,0.0E0)
                  IF (KODE.EQ.2) THEN
                    ZR=Z
                    CA=CSQRT(ZR)
                    ZR=ZR*CA*CMPLX(C23,0.0E0)
                    AA=REAL(ZR)
                    AA=ABS(AA)
                    CA = ZR - CMPLX(AA,0.0E0)
                    CV = CEXP(CA)*CV
                  ENDIF
                  CY = CW - CY - CV
                  ER(1) = CABS(CY)/CABS(CV)
                  JB = 1
                  JL = 1
                ELSE
C-----------------------------------------------------------------------
C     CHECKS IN -PI.LT.ARG(Z).LT.-PI/3 AND PI/3.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     CHECK AI    TEST #2
C-----------------------------------------------------------------------
                  CALL CAIRY(Z, 0, KODE, Y(2), NZ, IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  ZR=-Z
                  CV=CSQRT(ZR)
                  ZR=ZR*CV*CMPLX(C23,0.0E0)
                  CALL CBESJ(ZR,C23,KODE,2,YY,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CY=CMPLX(C43,0.0E0)*YY(1)/ZR - YY(2)
                  CA = YY(1)
                  CALL CBESJ(ZR,C13,KODE,2,YY,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  IF (KODE.EQ.2) THEN
                    AB = AIMAG(ZR)
                    AB = ABS(AB)
                    CW = CSQRT(Z)
                    ZW = Z*CW*CMPLX(C23,0.0E0)
                    CW = ZW+CMPLX(AB,0.0E0)
                    CW = CEXP(CW)
                    YY(1) = YY(1)*CW
                    YY(2) = YY(2)*CW
                    CY = CY*CW
                    CA = CA*CW
                    SC = CW
                  ENDIF
                  CW = CV*CMPLX(C13,0.0E0)
                  W(2) = CW*(YY(1)+CY)
                  ER(2) = CABS(Y(2)-W(2))
                  IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN
                    ER(2)=ER(2)/CABS(Y(2))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(2) = ER(2)/CABS(SC)
                    ENDIF
                  ENDIF
C-----------------------------------------------------------------------
C     CHECK AI'   TEST #3
C-----------------------------------------------------------------------
                  CY=CMPLX(C23,0.0E0)*YY(1)/ZR - YY(2)
                  W(3) = Z*CMPLX(C13,0.0E0)*(CY-CA)
                  CALL CAIRY(Z,1,KODE,Y(3),NZ,IERR)
                  ER(3) = CABS(Y(3)-W(3))
                  IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN
                    ER(3)=ER(3)/CABS(Y(3))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(3) = ER(3)/CABS(SC)
                    ENDIF
                  ENDIF
C-----------------------------------------------------------------------
C     CHECK BI    TEST #4
C-----------------------------------------------------------------------
                  CALL CBESH(ZR,C13,KODE,1,1,Y,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL CBESH(ZR,C13,KODE,2,1,YY,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CONA = CON1
                  CONB = CON2
                  CONC = CON2
                  COND = CON1
                  IF (KODE.EQ.2) THEN
                    AA = REAL(ZW)
                    AA = ABS(AA)
                    ZW =  CI*ZR - CMPLX(AA,0.0E0)
                    CW = CEXP(ZW)
                    CONA = CONA*CW
                    CONC = CONC*CW
                    ZW = -CI*ZR - CMPLX(AA,0.0E0)
                    CW = CEXP(ZW)
                    CONB = CONB*CW
                    COND = COND*CW
                    SC = CW
                  ENDIF
                  CW = CONA*Y(1)-CONB*YY(1)
                  CW = CV*CAV*CW
                  W(4) = CW*CHI
                  CALL CBIRY(Z,0,KODE,Y(4),IERR)
                  ER(4) = CABS(Y(4)-W(4))
                  IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN
                    ER(4)=ER(4)/CABS(Y(4))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(4) = ER(4)/CABS(SC)
                    ENDIF
                  ENDIF
C-----------------------------------------------------------------------
C     CHECK BI'   TEST #5
C-----------------------------------------------------------------------
                  CALL CBESH(ZR,C23,KODE,1,1,Y,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL CBESH(ZR,C23,KODE,2,1,YY,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CW = CONC*Y(1)-COND*YY(1)
                  CW = -Z*CAV*CW
                  W(5) = CW*CHI
                  CALL CBIRY(Z,1,KODE,Y(5),IERR)
                  ER(5) = CABS(Y(5)-W(5))
                  IF (ZY.NE.0.0D0.OR.ZX.GE.0.0D0) THEN
                    ER(5)=ER(5)/CABS(Y(5))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(5) = ER(5)/CABS(SC)
                    ENDIF
                  ENDIF
                  JB = 2
                  JL = 5
                ENDIF
                DO 190 J=JB,JL
                IF (ER(J).LT.ERTOL) GO TO 190
                IF (LFLG.EQ.1) GO TO 130
                write ( *,99995) ERTOL
99995           FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST W
     *ITH ERTOL =', E12.4/)
                write ( *,99994)
99994           FORMAT (/' OUTPUT FORMAT'/' KODE,IR,IT,IRSET,ICASE')
                write ( *,99993)
99993           FORMAT (' ER'/' I, Z, Y(I), W(I), ON THE I-TH TEST, I=1,
     *5'/)
                LFLG = 1
  130           CONTINUE
                write ( *,99992) KODE, IR, IT, IRSET, ICASE
99992           FORMAT (5I5)
                write ( *,99991) ER(J)
                write ( *,99989) J, Z, Y(J), W(J)
99991           FORMAT (E12.4)
99989           FORMAT (I5,6E12.4)
  190           CONTINUE
  140         CONTINUE
  150       CONTINUE
  160     CONTINUE
  170   CONTINUE
  180 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine ZQCBH

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C                *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCBH IS A QUICK CHECK ROUTINE FOR THE COMPLEX H BESSEL FUNCTIONS
C     GENERATED BY SUBROUTINE ZBESH.
C
C     ZQCBH GENERATES SEQUENCES OF H BESSEL FUNCTIONS FOR KIND 2 FROM
C     ZBESH AND CHECKS THEM AGAINST ANALYTIC CONTINUATION FORMULAS
C     IN THE (Z,FNU) SPACE:
C
C     KODE = 1 TESTS (ANALYTIC CONTINUATION FORMULAE, I**2 = -1):
C
C     H(FNU,2,Z)=-EXP(I*PI*FNU)*H(FNU,1,-Z),       -PI.LT.ARG(Z).LE.0
C
C               = 2*COS(PI*FNU)*H(FNU,2,-Z) + EXP(I*PI*FNU)*H(FNU,1,-Z),
C
C                                                   0.LT.ARG(Z).LE.PI
C
C     KODE = 2 TESTS FOR KINDS 1 AND 2:
C
C            EXP(-I*Z)*H(FNU,1,Z) = [EXP(-I*Z)*H(FNU,1,Z)]
C
C            EXP( I*Z)*H(FNU,2,Z) = [EXP( I*Z)*H(FNU,2,Z)]
C
C     WHERE THE LEFT SIDE IS COMPUTED WITH KODE = 1 AND THE RIGHT SIDE
C     WITH KODE = 2.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX CW, CI, U, V, W, Y, Z, ZN, CSGN
      EXTERNAL ZABS
      DOUBLE PRECISION AA, AB, AER, ALIM, ATOL, AV, CT, DIG, ERR, ELIM,
     * EPS, ER, ERTOL, FNU, FNUL, PI, R, RL, RM, D1M4, D1M5, R2, ST,
     * T, TOL, TS, XNU, D1MACH, SLAK, FILM, STR, STI, UR, UI, VR, VI,
     * WR, WI, YR, YI, CWR, CWI, CSGNR, CSGNI, ZR, ZI, ZNR, ZNI, ZABS
      INTEGER I, ICASE, IERR, IHP, IL, IR, IRB, IT, ITL, K, KODE, KK,
     *K1, K2, LFLG, LUN, MFLG, M, N, NU, NZ1, NZ2, NZ3, I1MACH, KEPS,
     *MQC, NL, NUL, KDO
      DIMENSION T(20), AER(20), XNU(20), UR(20), UI(20), VR(20), VI(20),
     * WR(20), WI(20), YR(20), YI(20), KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      D1M4 = D1MACH(4)
      TOL = MAX(D1M4,1.0D-18)
      AA = -DLOG10(D1M4)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      D1M5 = D1MACH(5)
      K = MIN(IABS(K1),IABS(K2))
      ELIM = 2.303D0*(DBLE(FLOAT(K))*D1M5-3.0D0)
      AB = AA*2.303D0
      ALIM = ELIM + MAX(-AB,-41.45D0)
      DIG = MIN(AA,18.0D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0
      SLAK = MAX(SLAK,3.0D0)
      ERTOL = TOL*10.0D0**SLAK
      RM = 0.5D0*(ALIM + ELIM)
      RM = MIN(RM,200.0D0)
      RM = MAX(RM,RL+10.0D0)
      R2 = MIN(FNUL,RM)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE H BESSEL FUNCTIONS FROM ZBES
     *H'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      ATOL = 100.0D0*TOL
      PI = 4.0D0*DATAN(1.0D0)
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0D0
        XNU(2) = 1.0D0
        XNU(3) = 2.0D0
        XNU(4) = 0.5D0*FNUL
        XNU(5) = FNUL + 1.1D0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0D0
        XNU(2) = 0.6D0
        XNU(3) = 1.3D0
        XNU(4) = 2.0D0
        XNU(5) = 0.5D0*FNUL
        XNU(6) = FNUL + 1.1D0
      ENDIF
      I = 2
      EPS = 0.01D0
      FILM=DBLE(FLOAT(IL-1))
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 170 KODE=1,2
        DO 160 N=1,NL
          DO 150 NU=1,NUL
            FNU = XNU(NU)
            DO 140 ICASE=1,3
              IRB = MIN(ICASE,2)
              DO 130 IR=IRB,3
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R =(EPS*DBLE(FLOAT(3-IR))+2.0D0*DBLE(FLOAT(IR-1)))/2.0D0
                GO TO 80
   60           CONTINUE
                R = (2.0D0*DBLE(FLOAT(3-IR))+R2*DBLE(FLOAT(IR-1)))/2.0D0
                GO TO 80
   70           CONTINUE
                IF (R2.GE.RM) GO TO 140
                R = (R2*DBLE(FLOAT(3-IR))+RM*DBLE(FLOAT(IR-1)))/2.0D0
   80           CONTINUE
                DO 120 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0D0
                  IF (ABS(ST).LT.ATOL) ST = 0.0D0
                  ZR = R*CT
                  ZI = R*ST
                  IF (KODE.EQ.1) THEN
                    M=2
                    CALL ZBESH(ZR,ZI,FNU,KODE,M,N,YR,YI,NZ1,IERR)
                    IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120
                    IF (ST.LT.0.0D0 .OR. (ST.EQ.0.0D0.AND.CT.GT.0.0D0))
     *              THEN
                      IHP = 1
                      ZNR = -ZR
                      ZNI = -ZI
                      M=1
                      CALL ZBESH(ZNR,ZNI,FNU,KODE,M,N,WR,WI,NZ2,IERR)
                      IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    ELSE
                      IHP = 2
                      ZNR = -ZR
                      ZNI = -ZI
                      M=2
                      CALL ZBESH(ZNR,ZNI,FNU,KODE,M,N,WR,WI,NZ3,IERR)
                      IF (IERR.NE.0.OR.NZ3.NE.0) GO TO 120
                      M=1
                      CALL ZBESH(ZNR,ZNI,FNU,KODE,M,N,VR,VI,NZ2,IERR)
                      IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    ENDIF
                    AB=MOD(FNU,2.0D0)*PI
                    CSGNR = COS(AB)
                    CSGNI = SIN(AB)
                    MFLG = 0
                    DO 100 I=1,N
                      AB = FNU+DBLE(FLOAT(I-1))
                      AA = MAX(0.5D0,AB)
                      IF (IHP.EQ.1) THEN
                        VR(I) = -(CSGNR*WR(I)-CSGNI*WI(I))
                        VI(I) = -(CSGNR*WI(I)+CSGNI*WR(I))
                        CWR = YR(I) - VR(I)
                        CWI = YI(I) - VI(I)
                      ELSE
                        CWR = CSGNR+CSGNR
                        STR =   CWR*WR(I) + CSGNR*VR(I)-CSGNI*VI(I)
                        VI(I) = CWR*WI(I) + CSGNR*VI(I)+CSGNI*VR(I)
                        VR(I) = STR
                        CWR = YR(I) - VR(I)
                        CWI = YI(I) - VI(I)
                      ENDIF
                      AV = ZABS(YR(I),YI(I))
                      ER = ZABS(CWR,CWI)
                      IF(ZI.EQ.0.0D0) THEN
                        IF(ABS(ZR).LT.AA) ER = ER/AV
                      ELSE
                        ER = ER/AV
                      ENDIF
                      AER(I) = ER
                      IF (ER.GT.ERTOL) MFLG = 1
                      CSGNR = -CSGNR
                      CSGNI = -CSGNI
  100               CONTINUE
                  ELSE
                    M=1
                    KK=1
                    CALL ZBESH(ZR,ZI,FNU,KK,M,N,UR,UI,NZ1,IERR)
                    IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120
                    CALL ZBESH(ZR,ZI,FNU,KODE,M,N,VR,VI,NZ2,IERR)
                    IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    M=2
                    KK=1
                    CALL ZBESH(ZR,ZI,FNU,KK,M,N,WR,WI,NZ1,IERR)
                    IF (IERR.NE.0.OR.NZ1.NE.0) GO TO 120
                    CALL ZBESH(ZR,ZI,FNU,KODE,M,N,YR,YI,NZ2,IERR)
                    IF (IERR.NE.0.OR.NZ2.NE.0) GO TO 120
                    ZNR = -ZI
                    ZNI =  ZR
                    CALL ZEXP(ZNR,ZNI,ZNR,ZNI)
                    MFLG = 0
                    DO 105 I=1,N
                      AB = FNU+DBLE(FLOAT(I-1))
                      AA = MAX(0.5D0,AB)
                      CALL ZDIV(UR(I),UI(I),ZNR,ZNI,STR,STI)
                      CWR = STR - VR(I)
                      CWI = STI - VI(I)
                      AV = ZABS(VR(I),VI(I))
                      ER = ZABS(CWR,CWI)
                      IF(ZI.EQ.0.0D0) THEN
                        IF(ABS(ZR).LT.AA) ER = ER/AV
                      ELSE
                        ER = ER/AV
                      ENDIF
                      ERR = ER
                      IF (ER.GT.ERTOL) MFLG = 1
                      CWR = ZNR*WR(I) - ZNI*WI(I) - YR(I)
                      CWI = ZNR*WI(I) + ZNI*WR(I) - YI(I)
                      AV = ZABS(YR(I),YI(I))
                      ER = ZABS(CWR,CWI)
                      IF(ZI.EQ.0.0D0) THEN
                        IF(ABS(ZR).LT.AA) ER = ER/AV
                      ELSE
                        ER = ER/AV
                      ENDIF
                      IF (ER.GT.ERTOL) MFLG = 1
                      AER(I) = ER+ERR
  105               CONTINUE
                  ENDIF
                  IF (MFLG.EQ.0) GO TO 120
                  IF (LFLG.EQ.1) GO TO 110
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', D12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,V(1),Y(1)')
                  LFLG = 1
  110             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, ICASE
99992             FORMAT (5I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) ZR,ZI,FNU,VR(1),VI(1),YR(1),YI(1)
99991             FORMAT (7D12.4)
  120           CONTINUE
  130         CONTINUE
  140       CONTINUE
  150     CONTINUE
  160   CONTINUE
  170 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine ZQCBI

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C                *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCBI IS A QUICK CHECK ROUTINE FOR THE COMPLEX I BESSEL FUNCTION
C     GENERATED BY SUBROUTINE ZBESI.
C
C     ZQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM
C     ZBESI AND CBESK AND CHECKS THE WRONSKIAN EVALUATION
C
C           I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z
C
C     IN THE RIGHT HALF PLANE AND A MODIFIED FORM
C
C          I(FNU+1,Z)*K(FNU,ZR) - I(FNU,Z)*K(FNU+1,ZR) = C/Z
C
C     IN THE LEFT HALF PLANE WHERE ZR=-Z AND C=EXP(I*FNU*SGN) WITH
C     SGN=+1 FOR IM(Z).GE.0 AND SGN=-1 FOR IM(Z).LT.0.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX  CONE, CSGN, CC, CV, CW, CY, W, Y, Z, ZR
      EXTERNAL ZABS
      DOUBLE PRECISION AA, AB, AER, ALIM, ARG, ATOL, CCR, CONER, CONEI,
     * CSGNI, CSGNR, CWI, CWR, CYI, CYR, CVR, CVI, DIG, ELIM, EPS, ER,
     * ERTOL, FNU, FNUL, GNU, HPI, PI, R, RL, R1M4, R1M5, R2, RM,
     * STI, STR, T, TOL, WI, WR, YI, YR, ZI, ZR, ZRI, ZRR, D1MACH, ZABS,
     * FILM, SLAK, TS, ST, CT, FFNU, XNU
      INTEGER I, ICASE, IL, IFNU, IPRNT, IR, IT, ITL, K, KK, KODE, K1,
     * K2, LFLG, LUN, MFLG, N, NZ, N1, NL, NU, NUL, MQC, IERR, I1MACH,
     * KEPS, KDO, IRB
      DIMENSION T(20), AER(20), YR(20), YI(20), WR(20), WI(20), XNU(20),
     *  KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = D1MACH(4)
      TOL = MAX(R1M4,1.0D-18)
      AA = -DLOG10(R1M4)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      R1M5 = D1MACH(5)
      K = MIN(IABS(K1),IABS(K2))
      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
      AB = AA*2.303D0
      ALIM = ELIM + MAX(-AB,-41.45D0)
      DIG = MIN(AA,18.0D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0
      SLAK = MAX(SLAK,3.0D0)
      ERTOL = TOL*10.0D0**SLAK
      RM = 0.5D0*(ALIM + ELIM)
      RM = MIN(RM,200.0D0)
      RM = MAX(RM,RL+10.0D0)
      R2 = MIN(RM,FNUL)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE I BESSEL FUNCTION FROM ZBESI
     *'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      CONER = 1.0D0
      CONEI = 0.0D0
      ATOL = 100.0D0*TOL
      HPI = 2.0D0*DATAN(1.0D0)
      PI = HPI + HPI
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0E0
        XNU(2) = 1.0E0
        XNU(3) = 2.0E0
        XNU(4) = 0.5E0*FNUL
        XNU(5) = FNUL + 1.1E0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0E0
        XNU(2) = 0.6E0
        XNU(3) = 1.3E0
        XNU(4) = 2.0E0
        XNU(5) = 0.5E0*FNUL
        XNU(6) = FNUL + 1.1E0
      ENDIF
      I = 2
      EPS = 0.01D0
      FILM=DBLE(FLOAT(IL-1))
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 200 KODE=1,2
        DO 190 N=1,NL
          N1 = N + 1
          DO 180 NU=1,NUL
            FNU = XNU(NU)
            IFNU = INT(SNGL(FNU))
            FFNU = FNU - DBLE(FLOAT(IFNU))
            ARG = PI*FFNU
            CSGNR = DCOS(ARG)
            CSGNI = DSIN(ARG)
            IF (MOD(IFNU,2).EQ.0) GO TO 50
            CSGNR = -CSGNR
            CSGNI = -CSGNI
   50       CONTINUE
            DO 170 ICASE=1,3
              IRB = MIN(2,ICASE)
              DO 160 IR=IRB,4
                GO TO (60, 70, 80), ICASE
   60           CONTINUE
                R = (0.2D0*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/
     *           3.0D0
                GO TO 90
   70           CONTINUE
                R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0
                GO TO 90
   80           CONTINUE
                IF (R2.EQ.RM) GO TO 180
                R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0
   90           CONTINUE
                DO 150 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (DABS(CT).LT.ATOL) CT = 0.0D0
                  IF (DABS(ST).LT.ATOL) ST = 0.0D0
                  ZR = R*CT
                  ZI = R*ST
                  IF (CT.GE.0.0E0) THEN
C-----------------------------------------------------------------------
C     WRONSKIAN CHECKS IN THE RIGHT HALF PLANE
C-----------------------------------------------------------------------
                    CALL ZBESI(ZR, ZI, FNU, KODE, N1, WR, WI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CALL ZBESK(ZR, ZI, FNU, KODE, N1, YR, YI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS
C     ON KODE=2
C-----------------------------------------------------------------------
                    CALL ZDIV(CONER,CONEI,ZR,ZI,CVR,CVI)
                    IF (KODE.EQ.2) THEN
                      STR = COS(ZI)
                      STI = SIN(ZI)
                      AA = STR*CVR - STI*CVI
                      CVI = STR*CVI + STI*CVR
                      CVR = AA
                    ENDIF
                    CCR = 1.0D0
                  ELSE
C-----------------------------------------------------------------------
C     WRONSKIAN CHECKS IN THE LEFT HALF PLANE
C-----------------------------------------------------------------------
                    ZRR = -ZR
                    ZRI = -ZI
                    CALL ZBESI(ZR, ZI, FNU, KODE, N1, YR, YI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CALL ZBESK(ZRR,ZRI, FNU, KODE, N1, WR, WI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CVR = CSGNR
                    CVI = CSGNI
                    IF (ZI.LT.0.0E0) THEN
                      CVI = -CVI
                    ENDIF
                    CALL ZDIV(CVR,CVI,ZR,ZI,CVR,CVI)
                    IF (KODE.EQ.2) THEN
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS
C     ON KODE=2. SCALE FACTOR = EXP(-I*YY) FOR RE(Z).LT.0
C-----------------------------------------------------------------------
                      CWR = COS(ZI)
                      CWI = -SIN(ZI)
                      STR = CVR*CWR-CVI*CWI
                      CVI = CVR*CWI+CVI*CWR
                      CVR = STR
                    ENDIF
                    CCR = -1.0D0
                  ENDIF
                  MFLG = 0
                  KK = 0
                  DO 130 I=1,N
                    CWR = WR(I)*YR(I+1)-WI(I)*YI(I+1)
                    CWI = WR(I)*YI(I+1)+WI(I)*YR(I+1)
                    CYR = YR(I)*WR(I+1)-YI(I)*WI(I+1)
                    CYI = YR(I)*WI(I+1)+YI(I)*WR(I+1)
                    CYR = CCR*CYR
                    CYI = CCR*CYI
                    CYR = CYR + CWR - CVR
                    CYI = CYI + CWI - CVI
                    ER = ZABS(CYR,CYI)/ZABS(CVR,CVI)
                    AER(I) = ER
                    IF (ER.GT.ERTOL) THEN
                      IF(KK.EQ.0) THEN
                        MFLG = 1
                        KK=I
                      ENDIF
                    ENDIF
                    IF (CT.LT.0.0E0) THEN
                      CVR = -CVR
                      CVI = -CVI
                    ENDIF
  130             CONTINUE
                  IF (MFLG.EQ.0) GO TO 150
                  IF (LFLG.EQ.1) GO TO 140
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK)        KK=INDEX O
     *F FIRST NON-ZERO PAIR'/)
                  LFLG = 1
  140             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, ICASE, KK
99992             FORMAT (6I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) ZR, ZI, FNU, YR(KK), YI(KK)
99991             FORMAT (6D12.4)
  150           CONTINUE
  160         CONTINUE
  170       CONTINUE
  180     CONTINUE
  190   CONTINUE
  200 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      IF (MQC.EQ.1) then
        return
      end if
C-----------------------------------------------------------------------
C     CHECKS NEAR UNDERFLOW LIMITS ON SERIES(I=1) AND UNIFORM
C     ASYMPTOTIC EXPANSION(I=2)
C-----------------------------------------------------------------------
      write ( *,99989)
99989 FORMAT (/' CHECKS NEAR UNDERFLOW AND OVERFLOW LIMITS'/)
      ZR = 1.4D0
      ZI = 1.4D0
      IPRNT = 0
      DO 280 I=1,2
        FNU = 10.2D0
        KODE = 1
        N = 20
  230   CONTINUE
        CALL ZBESI(ZR, ZI, FNU, KODE, N, YR, YI, NZ, IERR)
        IF (NZ.NE.0) GO TO 240
        FNU = FNU + 5.0D0
        GO TO 230
  240   CONTINUE
        IF (NZ.LT.10) GO TO 250
        FNU = FNU - 1.0D0
        GO TO 230
  250   CONTINUE
        CALL ZBESK(ZR, ZI, FNU, KODE, 2, WR, WI, NZ, IERR)
        CALL ZDIV(CONER, CONEI, ZR, ZI, STR, STI)
        CYR = WR(1)*YR(2) - WI(1)*YI(2)
        CYI = WR(1)*YI(2) + WI(1)*YR(2)
        CWR = WR(2)*YR(1) - WI(2)*YI(1)
        CWI = WR(2)*YI(1) + WI(2)*YR(1)
        CWR = CWR + CYR - STR
        CWI = CWI + CYI - STI
        ER = ZABS(CWR,CWI)/ZABS(STR,STI)
        IF (ER.LT.ERTOL) GO TO 270
        IF (IPRNT.EQ.1) GO TO 260
        write ( *,99988)
99988   FORMAT (/' OUTPUT FORMAT/19H ERROR,Z,FNU,KODE,N'/)
        IPRNT = 1
  260   CONTINUE
        write ( *,99987) ER, ZR, ZI, FNU, KODE, N
99987   FORMAT (4D12.4, 2I5)
  270   CONTINUE
        ZR = RL +RL
        ZI = 0.0D0
  280 CONTINUE
C-----------------------------------------------------------------------
C     CHECK NEAR OVERFLOW LIMITS
C-----------------------------------------------------------------------
      ZR = ELIM
      ZI = 0.0D0
      FNU = 0.0D0
  290 CONTINUE
      CALL ZBESK(ZR, ZI, FNU, KODE, N, YR, YI, NZ, IERR)
      IF (NZ.LT.10) GO TO 300
      IF (NZ.EQ.N) FNU = FNU + 3.0D0
      FNU = FNU + 2.0D0
      GO TO 290
  300 CONTINUE
      GNU = FNU + DBLE(FLOAT(N-2))
      CALL ZBESI(ZR, ZI, GNU, KODE, 2, WR, WI, NZ, IERR)
      CALL ZDIV(CONER, CONEI, ZR, ZI, STR, STI)
      CYR = YR(N-1)*WR(2) - YI(N-1)*WI(2)
      CYI = YR(N-1)*WI(2) + YI(N-1)*WR(2)
      CWR = YR(N)*WR(1) - YI(N)*WI(1)
      CWI = YR(N)*WI(1) + YI(N)*WR(1)
      CWR = CWR + CYR - STR
      CWI = CWI + CYI - STI
      ER = ZABS(CWR,CWI)/ZABS(STR,STI)
      IF (ER.LT.ERTOL) GO TO 320
      IF (IPRNT.EQ.1) GO TO 310
      write ( *,99988)
      IPRNT = 1
  310 CONTINUE
      write ( *,99987) ER, ZR, ZI, FNU, KODE, N
  320 CONTINUE
      IF (IPRNT.EQ.0) write ( *,99990)
      return
      END
      subroutine ZQCBJ

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C                *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCBJ IS A QUICK CHECK ROUTINE FOR THE COMPLEX J BESSEL FUNCTION
C     GENERATED BY SUBROUTINE ZBESJ.
C
C     ZQCBJ GENERATES SEQUENCES OF J AND H BESSEL FUNCTIONS FROM ZBESJ
C     AND ZBESH AND CHECKS THE WRONSKIANS
C
C     J(FNU,Z)*H(FNU+1,1,Z)-J(FNU+1,Z)*H(FNU,1,Z)=2/(PI*I*Z)   Y.GE.0
C
C     J(FNU,Z)*H(FNU+1,2,Z)-J(FNU+1,Z)*H(FNU,2,Z)=-2/(PI*I*Z)  Y.LT.0
C
C     IN THEIR RESPECTIVE HALF PLANES.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX Z, WR, CJ, CH, CON, T1, T2, CER
      EXTERNAL ZABS
      DOUBLE PRECISION AA, AB, AER, ALIM, ATOL, CONR, CONI, WRR, WRI,
     * CERR, CERI, CT, DIG, ELIM, EPS, ER, ERTOL, FNU, FNUL, GNU, HPI,
     * PI, R, RL, RM, R1M4, R1M5, R2, ST, STR, STI, TOL, T1R, T1I, T2R,
     * T2I, T, XNU, ZI, ZR, D1MACH, ZABS, FILM, SLAK, TS, SGN, CJR, CJI,
     * CHR, CHI
      INTEGER I, ICASE, IL, IR, IRB, IT, ITL, K, KK, KODE, K1, K2,
     * LFLG, LUN, M, MFLG, N, NU, NZJ, NZH, I1MACH, KEPS, KDO, NL,
     *NUL, MQC, IERRJ, IERRH
      DIMENSION T(20), AER(25), XNU(20), CJR(20), CJI(20), CHR(20),
     * CHI(20), KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = D1MACH(4)
      TOL = DMAX1(R1M4,1.0D-18)
      AA = -DLOG10(R1M4)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      R1M5 = D1MACH(5)
      K = MIN0(ABS(K1),ABS(K2))
      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
      AB = AA*2.303D0
      ALIM = ELIM + DMAX1(-AB,-41.45D0)
      DIG = DMIN1(AA,18.0D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0
      SLAK = DMAX1(SLAK,3.0D0)
      ERTOL = TOL*10.0D0**SLAK
      RM = 0.5D0*(ALIM + ELIM)
      RM = DMIN1(RM,200.0D0)
      RM = DMAX1(RM,RL+10.0D0)
      R2 = DMIN1(RM,FNUL)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE J BESSEL FUNCTION FROM ZBESJ
     *'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      ATOL = 100.0D0*TOL
      HPI = 2.0D0*DATAN(1.0D0)
      PI = HPI + HPI
      CONR = 0.0D0
      CONI = -1.0D0/HPI
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0D0
        XNU(2) = 1.0D0
        XNU(3) = 2.0D0
        XNU(4) = 0.5D0*FNUL
        XNU(5) = FNUL + 1.1D0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0D0
        XNU(2) = 0.6D0
        XNU(3) = 1.3D0
        XNU(4) = 2.0D0
        XNU(5) = 0.5D0*FNUL
        XNU(6) = FNUL + 1.1D0
      ENDIF
      I = 2
      EPS = 0.01D0
      FILM=DBLE(FLOAT(IL-1))
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 260 KODE=1,2
        DO 250 N=1,NL
          NP = N+1
          DO 240 NU=1,NUL
            FNU = XNU(NU)
            GNU = FNU + DBLE(FLOAT(N-1)) + 1.0D0
            GNU = SQRT(GNU)
            GNU = MIN(GNU, 0.5D0*RL)
            DO 230 ICASE=1,3
              IRB = MIN0(2,ICASE)
              DO 220 IR=IRB,4
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (GNU*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/
     *           3.0D0
                GO TO 80
   60           CONTINUE
                R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0
                GO TO 80
   70           CONTINUE
                IF (R2.GE.RM) GO TO 230
                R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0
   80           CONTINUE
                DO 210 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (ABS(CT).LT.ATOL) CT = 0.0D0
                  IF (ABS(ST).LT.ATOL) ST = 0.0D0
                  ZR = R*CT
                  ZI = R*ST
                  IF(ZR.EQ.0.0.AND.ZI.EQ.0.0) GO TO 210
                  CALL ZDIV(CONR,CONI,ZR,ZI,WRR,WRI)
                  M=1
                  IF(ZI.LT.0.0) THEN
                    M=2
                    WRR = -WRR
                    WRI = -WRI
                  ENDIF
                  CALL ZBESJ(ZR,ZI,FNU,KODE,NP,CJR,CJI,NZJ,IERRJ)
                  CALL ZBESH(ZR,ZI,FNU,KODE,M,NP,CHR,CHI,NZH,IERRH)
                  IF(NZJ.NE.0.OR.NZH.NE.0) GO TO 210
                  IF(IERRJ.NE.0.OR.IERRH.NE.0) GO TO 210
                  IF(KODE.EQ.2) THEN
                    SGN = 3.0D0-2.0D0*DBLE(FLOAT(M))
                    STR = COS(ZR)
                    STI = -SGN*SIN(ZR)
                    T1R = WRR*STR - WRI*STI
                    WRI = WRR*STI + WRI*STR
                    WRR = T1R
                  ENDIF
                  KK = 0
                  MFLG = 0
                  DO 190 I=1,N
                    STR = CJR(I)*CHR(I+1) - CJI(I)*CHI(I+1)
                    T1I = CJR(I)*CHI(I+1) + CJI(I)*CHR(I+1)
                    T1R = STR
                    STR = CJR(I+1)*CHR(I) - CJI(I+1)*CHI(I)
                    T2I = CJR(I+1)*CHI(I) + CJI(I+1)*CHR(I)
                    T2R = STR
                    CERR = T1R - T2R -WRR
                    CERI = T1I - T2I -WRI
                    ER=ZABS(CERR,CERI)/ZABS(WRR,WRI)
                    IF (ER.GT.ERTOL) THEN
                      IF(MFLG.EQ.0) THEN
                        MFLG = 1
                        KK=I
                      ENDIF
                    ENDIF
                    AER(I)=ER
  190             CONTINUE
                  IF (MFLG.EQ.0) GO TO 210
                  IF (LFLG.EQ.1) GO TO 200
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', D12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZJ,NZH,ICASE
     *')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,CJ(KK),CH(KK), KK=INDEX
     * OF FIRST NON-ZERO Y,W PAIR'/)
                  LFLG = 1
  200             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, NZJ, NZH, ICASE
99992             FORMAT (8I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) ZR, ZI, FNU, CJR(KK), CJI(KK),
     *             CHR(KK), CHI(KK)
99991             FORMAT (9D12.4)
  210           CONTINUE
  220         CONTINUE
  230       CONTINUE
  240     CONTINUE
  250   CONTINUE
  260 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine ZQCBK

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C                *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCBK IS A QUICK CHECK ROUTINE FOR THE COMPLEX K BESSEL FUNCTION
C     GENERATED BY SUBROUTINE ZBESK.
C
C     ZQCBK GENERATES SEQUENCES OF I AND K BESSEL FUNCTIONS FROM
C     ZBESI AND ZBESK AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION
C
C           I(FNU,Z)*K(FNU+1,Z) + I(FNU+1,Z)*K(FNU,Z) = 1/Z
C
C     IN THE RIGHT HALF PLANE AND THE ANALYTIC CONTINUATION FORMULA
C     FOR H(FNU,2,Z) IN TERMS OF THE K FUNCTION
C
C           K(FNU,Z) = C3*H(FNU,2,ZR) + C4*H(FNU,1,ZR)    IM(Z).GE.0
C
C                    = CONJG(K(FNU,CONJG(Z)))             IM(Z).LT.0
C
C     IN THE LEFT HALF PLANE WHERE C3=C1*CONJG(C2)*C5, C4 = C2*C5
C     C1=2*COS(PI*FNU),   C2=EXP(PI*FNU*I/2),   C5 =-PI*I/2   AND
C     ZR = Z*EXP(-3*PI*I/2) = Z*I
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX CONE, CSGN, CV, CW, CY, C1, C2, C3, C4, V, W, Y, Z, ZR,
C    * CIP
      EXTERNAL ZABS
      DOUBLE PRECISION AA, AB, AER, ALIM, ARG, ATOL, CONEI, CONER,
     * CSGNI, CSGNR, CVI, CVR, CWI, CWR, CYI, CYR, DIG, ELIM, EPS,
     * ER, ERTOL, FFNU, FNU, FNUL, HPI, PI, R, RL, RM, R1M4, R1M5, R2,
     * STI, STR, T, TOL, WI, WR, XNU, YI, YR, ZI, C1R, C1I, C2R, C2I,
     * ZRI, ZRR, ZR, D1MACH, ZABS, FILM, CT, ST, TS, SLAK, C3R, C3I,
     * C4R, C4I, VR, VI, CIPR, CIPI, COEI, ZZR, ZZI
      INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, I4, K, KK, KODE,
     * K1, K2, LFLG, LUN, M, MFLG, N, NU, N1, IERR, I1MACH,
     *KEPS, KDO, NL, NUL, NZ, MQC
      DIMENSION T(20), AER(25), XNU(20), VR(20), VI(20), YR(20), YI(20),
     * WR(20), WI(20), KEPS(20), KDO(20), CIPR(4), CIPI(4)

      DATA CIPR(1), CIPI(1), CIPR(2), CIPI(2), CIPR(3),CIPI(3),CIPR(4),
     * CIPI(4) / 1.0D0,0.0D0,0.0D0,1.0D0,-1.0D0,0.0D0,0.0D0,-1.0D0/
      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = D1MACH(4)
      TOL = MAX(R1M4,1.0D-18)
      AA = -DLOG10(R1M4)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      R1M5 = D1MACH(5)
      K = MIN(IABS(K1),IABS(K2))
      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
      AB = AA*2.303D0
      ALIM = ELIM + MAX(-AB,-41.45D0)
      DIG = MIN(AA,18.0D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0
      SLAK = MAX(SLAK,3.0D0)
      ERTOL = TOL*10.0D0**SLAK
      RM = 0.5D0*(ALIM + ELIM)
      RM = MIN(RM,200.0D0)
      RM = MAX(RM,RL+10.0D0)
      R2 = MIN(RM,FNUL)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (54H QUICK CHECK ROUTINE FOR THE K BESSEL FUNCTION FROM ZB,
     * 3HESK/)
      write ( *,99998)
99998 FORMAT (37H PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG)
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      CONER = 1.0D0
      CONEI = 0.0D0
      ATOL = 100.0D0*TOL
      HPI = 2.0D0*DATAN(1.0D0)
      PI = HPI + HPI
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI NEAR FORMULA BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        NUL=5
        XNU(1) = 0.0D0
        XNU(2) = 1.0D0
        XNU(3) = 2.0D0
        XNU(4) = 0.5D0*FNUL
        XNU(5) = FNUL + 1.1D0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(12)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        KEPS(10)=1
        KEPS(11)=1
        NUL=6
        XNU(1) = 0.0D0
        XNU(2) = 0.6D0
        XNU(3) = 1.3D0
        XNU(4) = 2.0D0
        XNU(5) = 0.5D0*FNUL
        XNU(6) = FNUL + 1.1D0
      ENDIF
      I = 2
      EPS = 0.01D0
      FILM=DBLE(FLOAT(IL-1))
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 200 KODE=1,2
        DO 190 N=1,NL
          N1 = N + 1
          DO 180 NU=1,NUL
            FNU = XNU(NU)
            IFNU = INT(SNGL(FNU))
            FFNU = FNU - DBLE(FLOAT(IFNU))
            ARG = HPI*FFNU
            CSGNR = DCOS(ARG)
            CSGNI = DSIN(ARG)
            I4 = MOD(IFNU,4)+1
            STR = CSGNR*CIPR(I4) - CSGNI*CIPI(I4)
            CSGNI = CSGNR*CIPI(I4) + CSGNI*CIPR(I4)
            CSGNR = STR
   50       CONTINUE
            DO 170 ICASE=1,3
              IRB = MIN(2,ICASE)
              DO 160 IR=IRB,4
                GO TO (60, 70, 80), ICASE
   60           CONTINUE
                R = (0.2D0*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/
     *           3.0D0
                GO TO 90
   70           CONTINUE
                R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0
                GO TO 90
   80           CONTINUE
                IF (R2.EQ.RM) GO TO 180
                R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0
   90           CONTINUE
                DO 150 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (DABS(CT).LT.ATOL) CT = 0.0D0
                  IF (DABS(ST).LT.ATOL) ST = 0.0D0
                  ZR = R*CT
                  ZI = R*ST
                  IF (ZR.GE.0.0E0) THEN
C-----------------------------------------------------------------------
C     WRONSKIAN CHECKS IN THE RIGHT HALF PLANE
C-----------------------------------------------------------------------
                    CALL ZBESI(ZR, ZI, FNU, KODE, N1, WR, WI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    CALL ZBESK(ZR, ZI, FNU, KODE, N1, YR, YI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO WRONSKIAN DUE TO SCALING OF I AND K FUNCTIONS
C     ON KODE=2
C-----------------------------------------------------------------------
                    CALL ZDIV(CONER,CONEI,ZR,ZI,CVR,CVI)
                    IF (KODE.EQ.2) THEN
                      STR = COS(ZI)
                      STI = SIN(ZI)
                      AA = STR*CVR - STI*CVI
                      CVI = STR*CVI + STI*CVR
                      CVR = AA
                    ENDIF
                    MFLG = 0
                    KK = 0
                    DO 130 I=1,N
                      CWR = WR(I)*YR(I+1)-WI(I)*YI(I+1)
                      CWI = WR(I)*YI(I+1)+WI(I)*YR(I+1)
                      CYR = YR(I)*WR(I+1)-YI(I)*WI(I+1)
                      CYI = YR(I)*WI(I+1)+YI(I)*WR(I+1)
                      CYR = CYR + CWR - CVR
                      CYI = CYI + CWI - CVI
                      ER = ZABS(CYR,CYI)/ZABS(CVR,CVI)
                      AER(I) = ER
                      IF (ER.GT.ERTOL) THEN
                        IF(KK.EQ.0) THEN
                          MFLG = 1
                          KK=I
                        ENDIF
                      ENDIF
  130               CONTINUE
                  ELSE
C-----------------------------------------------------------------------
C     ANALYTIC CONTINUATION FORMULA CHECKS FOR LEFT HALF PLANE IN TERMS
C     OF H(FNU,1,Z) AND H(FNU,2,Z)
C-----------------------------------------------------------------------
                    ZZR = ZR
                    ZZI = ZI
                    IF (ZI.LT.0.0E0) THEN
                      ZZI = -ZZI
                    ENDIF
                    ZRR = -ZZI
                    ZRI = ZZR
                    M=1
                    CALL ZBESH(ZRR,ZRI,FNU,KODE,M,N,WR,WI,NZ,IERR)
                    IF (IERR.NE.0) GO TO 150
                    M=2
                    CALL ZBESH(ZRR,ZRI,FNU,KODE,M,N,VR,VI,NZ,IERR)
                    IF (IERR.NE.0) GO TO 150
                    CALL ZBESK(ZR, ZI, FNU, KODE, N, YR, YI, NZ, IERR)
                    IF (NZ.NE.0.OR.IERR.NE.0) GO TO 150
                    COEI = -HPI
                    MFLG = 0
                    KK = 0
                    AA = 2.0D0*COS(PI*FFNU)
                    IF(MOD(IFNU,2).NE.0) AA = -AA
                    C1R = AA
                    C1I = 0.0D0
                    C2R = CSGNR
                    C2I = CSGNI
                    DO 135 I=1,N
                      C3R = C1R
                      C3I = C1I
                      C4R = C2R
                      C4I = C2I
                      IF (KODE.EQ.2) THEN
C-----------------------------------------------------------------------
C     ADJUSTMENTS TO COEFICIENTS DUE TO SCALING OF H(FNU,1,Z) AND
C     H(FNU,2,Z) FUNCTIONS ON KODE = 2.
C-----------------------------------------------------------------------
                        AB=ZABS(VR(I),VI(I))
                        AA=LOG(AB)+ZR+ZR
                        IF (AA.GT.ELIM) GO TO 150
                        IF (AA.LT.-ELIM) THEN
                          C3R = 0.0D0
                          C3I = 0.0D0
                        ELSE
                          STR = ZZR+ZZR
                          STI = ZZI+ZZI
                          CALL ZEXP(STR,STI,CVR,CVI)
                          STR = C3R*CVR - C3I*CVI
                          C3I = C3R*CVI + C3I*CVR
                          C3R = STR
                        ENDIF
                      ENDIF
C                     CY = (C3*CONJG(C2)*V(I)+C4*W(I))*COE
                      STR = C3R*C2R + C3I*C2I
                      STI = -C3R*C2I + C3I*C2R
                      CYR = STR*VR(I) - STI*VI(I)
                      CYI = STR*VI(I) + STI*VR(I)
                      CYR = CYR + (C4R*WR(I) - C4I*WI(I))
                      CYI = CYI + (C4R*WI(I) + C4I*WR(I))
                      STR = -CYI*COEI
                      CYI =  CYR*COEI
                      CYR = STR
                      IF (ZI.LT.0.0D0) THEN
                        CYI = -CYI
                      ENDIF
                      STR = CYR - YR(I)
                      STI = CYI - YI(I)
                      ER = ZABS(STR,STI)/ZABS(YR(I),YI(I))
                      AER(I) = ER
                      IF (ER.GT.ERTOL) THEN
                        IF(KK.EQ.0) THEN
                          MFLG = 1
                          KK=I
                        ENDIF
                      ENDIF
                      STR = -C2I
                      C2I = C2R
                      C2R = STR
                      C1R = -C1R
                      C1I = -C1I
  135               CONTINUE
                  ENDIF
                  IF (MFLG.EQ.0) GO TO 150
                  IF (LFLG.EQ.1) GO TO 140
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL =', E12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,ICASE,KK')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,Y(KK)        KK=INDEX O
     *F FIRST NON-ZERO PAIR'/)
                  LFLG = 1
  140             CONTINUE
                  write ( *,99992) KODE, N, IR, IT, ICASE, KK
99992             FORMAT (6I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) ZR, ZI, FNU, YR(KK), YI(KK)
99991             FORMAT (6D12.4)
  150           CONTINUE
  160         CONTINUE
  170       CONTINUE
  180     CONTINUE
  190   CONTINUE
  200 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      subroutine ZQCBY

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C                *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCBY IS A QUICK CHECK ROUTINE FOR THE COMPLEX Y BESSEL FUNCTION
C     GENERATED BY SUBROUTINE ZBESY.
C
C     ZQCBY GENERATES SEQUENCES OF Y BESSEL FUNCTIONS FROM ZBESY AND
C     ZBESYH AND COMPARES THEM FOR A VARIETY OF VALUES IN THE (Z,FNU)
C     SPACE. ZBESYH IS AN OLD VERSION OF ZBESY WHICH COMPUTES THE Y
C     FUNCTION FROM THE H FUNCTIONS OF KINDS 1 AND 2.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX CW,CWRK,V,W,Z
      EXTERNAL ZABS
      DOUBLE PRECISION AA, AB, AER, ALIM, ATOL, AV, CWI, CWR, CWRKI,
     * CWRKR, DIG, ELIM, EPS, ER, ERTOL, FFNU, FNU, FNUL, HPI, PI, R,
     * RL, RM, R1M4, R1M5, R2, T, TOL, VI, VR, WI, WR, XNU, ZI, ZR,
     * D1MACH, ZABS, ST, CT, TS, SLAK, FILM
      INTEGER I, ICASE, IFNU, IL, IR, IRB, IT, ITL, K, KK,
     * KODE, K1, K2, LFLG, LUN, MFLG, N, NU, NZ1, NZ2, IERR, I1MACH,
     * MQC, KDO, KEPS, NL, NUL
      DIMENSION T(20), AER(20), XNU(20), WR(20), WI(20), VR(20), VI(20),
     * CWRKR(20), CWRKI(20), KDO(20), KEPS(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = D1MACH(4)
      TOL = MAX(R1M4,1.0D-18)
      AA = -DLOG10(R1M4)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      R1M5 = D1MACH(5)
      K = MIN(IABS(K1),IABS(K2))
      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
      AB = AA*2.303D0
      ALIM = ELIM + MAX(-AB,-41.45D0)
      DIG = MIN(AA,18.0D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0
      SLAK = MAX(SLAK,3.0D0)
      ERTOL = TOL*10.0D0**SLAK
      RM = 0.5D0*(ALIM + ELIM)
      RM = MIN(RM,200.0D0)
      RM = MAX(RM,RL+10.0D0)
      R2 = MIN(RM,FNUL)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM ZBESY
     *'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      ATOL = 100.0D0*TOL
      HPI = 2.0D0*DATAN(1.0D0)
      PI = HPI + HPI
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI/2.LT.ARG(Z).LE.PI NEAR BOUNDARIES
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        NL=2
        IL=5
        DO 5 I=1,IL
          KEPS(I)=0
          KDO(I)=0
    5   CONTINUE
        KDO(5)=1
        NUL=5
        XNU(1) = 0.0D0
        XNU(2) = 1.0D0
        XNU(3) = 2.0D0
        XNU(4) = 0.5D0*FNUL
        XNU(5) = FNUL + 1.2D0
      ELSE
        NL=4
        IL=13
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KDO(2)=1
        KDO(6)=1
        KDO(8)=1
        KDO(11)=1
        KDO(12)=1
        KDO(13)=1
        KEPS(3)=1
        KEPS(4)=1
        KEPS(5)=1
        KEPS(9)=1
        NUL=6
        XNU(1) = 0.0D0
        XNU(2) = 0.6D0
        XNU(3) = 1.3D0
        XNU(4) = 2.0D0
        XNU(5) = 0.5D0*FNUL
        XNU(6) = FNUL + 1.2D0
      ENDIF
      I = 2
      EPS = 0.01D0
      FILM=DBLE(FLOAT(IL-1))
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF (KDO(K).EQ.0) THEN
          T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 190 KODE=1,2
        DO 180 N=1,NL
          DO 170 NU=1,NUL
            FNU = XNU(NU)
            IFNU = INT(SNGL(FNU))
            FFNU = FNU - DBLE(FLOAT(IFNU))
            DO 160 ICASE=1,3
              IRB = MIN(2,ICASE)
              DO 150 IR=IRB,4
                GO TO (50, 60, 70), ICASE
   50           CONTINUE
                R = (EPS*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/
     *           3.0D0
                GO TO 80
   60           CONTINUE
                R = (2.0D0*DBLE(FLOAT(4-IR))+R2*DBLE(FLOAT(IR-1)))/3.0D0
                GO TO 80
   70           CONTINUE
                IF (RM.EQ.R2) GO TO 160
                R = (R2*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0
   80           CONTINUE
                DO 140 IT=1,ITL
                  CT = COS(T(IT))
                  ST = SIN(T(IT))
                  IF (DABS(CT).LT.ATOL) CT = 0.0D0
                  IF (DABS(ST).LT.ATOL) ST = 0.0D0
                  ZR = R*CT
                  ZI = R*ST
                  CALL ZBESY(ZR, ZI, FNU, KODE, N, VR, VI, NZ2, CWRKR,
     &                       CWRKI, IERR)
                  IF (NZ2.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL ZBESYH(ZR, ZI, FNU, KODE, N, WR, WI, NZ1, CWRKR,
     &                        CWRKI, IERR)
                  IF (NZ1.NE.0.OR.IERR.NE.0) GO TO 140
                  MFLG = 0
                  DO 120 I=1,N
                    AB = FNU+DBLE(FLOAT(I-1))
                    AA = MAX(0.5D0,AB)
                    CWR = WR(I) - VR(I)
                    CWI = WI(I) - VI(I)
                    AV = ZABS(VR(I),VI(I))
                    ER = ZABS(CWR,CWI)
                    IF (AV.NE.0.0D0) THEN
                      IF (ZI.EQ.0.0D0) THEN
                        IF (ZR.GT.0.0D0) THEN
                          IF (DABS(ZR).LT.AA) ER = ER/AV
                        ELSE
                          IF (DABS(FFNU-0.5D0).LT.0.125D0) THEN
                            IF (DABS(ZR).LT.AA) ER = ER/AV
                          ELSE
                            ER = ER/AV
                          ENDIF
                        ENDIF
                      ELSE
                        ER = ER/AV
                      ENDIF
                    ENDIF
                    AER(I) = ER
                    IF (ER.GT.ERTOL) MFLG = 1
  120             CONTINUE
                  IF (MFLG.EQ.0) GO TO 140
                  IF (LFLG.EQ.1) GO TO 130
                  write ( *,99995) ERTOL
99995             FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST
     * WITH ERTOL = ', D12.4/)
                  write ( *,99994)
99994             FORMAT (/' OUTPUT FORMAT'/' KODE,N,IR,IT,NZ1,NZ2,ICASE
     *')
                  write ( *,99993)
99993             FORMAT (' ER(K),K=1,N'/' Z,FNU,W(KK),V(KK), KK=INDEX O
     *F FIRST NON-ZERO W,V PAIR'/)
                  LFLG = 1
  130             CONTINUE
                  KK = MAX(NZ1,NZ2) + 1
                  KK = MIN(N,KK)
                  write ( *,99992) KODE, N, IR, IT, NZ1, NZ2, ICASE
99992             FORMAT (8I5)
                  write ( *,99991) (AER(K),K=1,N)
                  write ( *,99991) ZR, ZI, FNU, WR(KK), WI(KK),
     *             VR(KK), VI(KK)
99991             FORMAT (7D12.4)
  140           CONTINUE
  150         CONTINUE
  160       CONTINUE
  170     CONTINUE
  180   CONTINUE
  190 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
      SUBROUTINE ZBESYH(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, CWRKR,
     *           CWRKI, IERR)

c*********************************************************************72
c
C***BEGIN PROLOGUE  ZBESYH
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C***CATEGORY NO.  B5K
C***KEYWORDS  Y-BESSEL FUNCTION,BESSEL FUNCTION OF COMPLEX ARGUMENT,
C             BESSEL FUNCTION OF SECOND KIND
C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
C***PURPOSE  TO COMPUTE THE Y-BESSEL FUNCTION OF A COMPLEX ARGUMENT
C***DESCRIPTION
C
C                ***A DOUBLE PRECISION ROUTINE***
C
C         ON KODE=1, ZBESYH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
C         BESSEL FUNCTIONS CY(I)=Y(FNU+I-1,Z) FOR REAL, NONNEGATIVE
C         ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
C         -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESYH RETURNS THE SCALED
C         FUNCTIONS
C
C         CY(I)=EXP(-ABS(Y))*Y(FNU+I-1,Z)   I = 1,...,N , Y=AIMAG(Z)
C
C         WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
C         LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
C         ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
C         (REF. 1).
C
C         INPUT      ZR,ZI,FNU ARE DOUBLE PRECISION
C           ZR,ZI  - Z=CMPLX(ZR,ZI), Z.NE.CMPLX(0.0D0,0.0D0),
C                    -PI.LT.ARG(Z).LE.PI
C           FNU    - ORDER OF INITIAL Y FUNCTION, FNU.GE.0.0D0
C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
C                    KODE= 1  RETURNS
C                             CY(I)=Y(FNU+I-1,Z), I=1,...,N
C                        = 2  RETURNS
C                             CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...,N
C                             WHERE Y=AIMAG(Z)
C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
C           CWRKR, - DOUBLE PRECISION WORK VECTORS OF DIMENSION AT
C           CWRKI    AT LEAST N
C
C         OUTPUT     CYR,CYI ARE DOUBLE PRECISION
C           CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
C                    CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
C                    CY(I)=Y(FNU+I-1,Z)  OR
C                    CY(I)=Y(FNU+I-1,Z)*EXP(-ABS(Y))  I=1,...,N
C                    DEPENDING ON KODE.
C           NZ     - NZ=0 , A NORMAL RETURN
C                    NZ.GT.0 , NZ COMPONENTS OF CY SET TO ZERO DUE TO
C                    UNDERFLOW (GENERALLY ON KODE=2)
C           IERR   - ERROR FLAG
C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
C                    IERR=1, INPUT ERROR   - NO COMPUTATION
C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU IS
C                            TOO LARGE OR CABS(Z) IS TOO SMALL OR BOTH
C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
C                            ACCURACY
C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
C                            CANCE BY ARGUMENT REDUCTION
C                    IERR=5, ERROR              - NO COMPUTATION,
C                            ALGORITHM TERMINATION CONDITION NOT MET
C
C***LONG DESCRIPTION
C
C         THE COMPUTATION IS CARRIED OUT BY THE FORMULA
C
C              Y(FNU,Z)=0.5*(H(1,FNU,Z)-H(2,FNU,Z))/I
C
C         WHERE I**2 = -1 AND THE HANKEL BESSEL FUNCTIONS H(1,FNU,Z)
C         AND H(2,FNU,Z) ARE CALCULATED IN ZBESH.
C
C         FOR NEGATIVE ORDERS,THE FORMULA
C
C              Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
C
C         CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO HALF ODD
C         INTEGERS THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE
C         POSITIVE HALF ODD INTEGER,THE MAGNITUDE OF Y(-FNU,Z)=J(FNU,Z)*
C         SIN(PI*FNU) IS A LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS
C         NOT A HALF ODD INTEGER, Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A
C         LARGE POSITIVE POWER OF TEN AND THE MOST THAT THE SECOND TERM
C         CAN BE REDUCED IS BY UNIT ROUNDOFF FROM THE COEFFICIENT. THUS,
C         WIDE CHANGES CAN OCCUR WITHIN UNIT ROUNDOFF OF A LARGE HALF
C         ODD INTEGER. HERE, LARGE MEANS FNU.GT.CABS(Z).
C
C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
C         IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
C         DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
C
C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
C         OR -PI/2+P.
C
C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
C                 COMMERCE, 1955.
C
C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
C
C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
C
C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
C                 1018, MAY, 1985
C
C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
C                 MATH. SOFTWARE, 12, NO. 3, SEPTEMBER 1986, PP 265-273.
C
C***ROUTINES CALLED  ZBESH,I1MACH,D1MACH
C***END PROLOGUE  ZBESYH
C
C     COMPLEX CWRK,CY,C1,C2,EX,HCI,Z,ZU,ZV
      DOUBLE PRECISION CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2R,
     * ELIM, EXI, EXR, EY, FNU, HCII, STI, STR, TAY, ZI, ZR, DEXP,
     * D1MACH, ASCLE, RTOL, ATOL, AA, BB, TOL, R1M5
      INTEGER I, IERR, K, KODE, K1, K2, N, NZ, NZ1, NZ2, I1MACH
      DIMENSION CYR(N), CYI(N), CWRKR(N), CWRKI(N)
C***FIRST EXECUTABLE STATEMENT  ZBESYH
      IERR = 0
      NZ=0
      IF (ZR.EQ.0.0D0 .AND. ZI.EQ.0.0D0) IERR=1
      IF (FNU.LT.0.0D0) IERR=1
      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
      IF (N.LT.1) IERR=1
      IF (IERR.NE.0) RETURN
      HCII = 0.5D0
      CALL ZBESH(ZR, ZI, FNU, KODE, 1, N, CYR, CYI, NZ1, IERR)
      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
      CALL ZBESH(ZR, ZI, FNU, KODE, 2, N, CWRKR, CWRKI, NZ2, IERR)
      IF (IERR.NE.0.AND.IERR.NE.3) GO TO 170
      NZ = MIN0(NZ1,NZ2)
      IF (KODE.EQ.2) GO TO 60
      DO 50 I=1,N
        STR = CWRKR(I) - CYR(I)
        STI = CWRKI(I) - CYI(I)
        CYR(I) = -STI*HCII
        CYI(I) = STR*HCII
   50 CONTINUE
      RETURN
   60 CONTINUE
      TOL = DMAX1(D1MACH(4),1.0D-18)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      K = MIN0(IABS(K1),IABS(K2))
      R1M5 = D1MACH(5)
C-----------------------------------------------------------------------
C     ELIM IS THE APPROXIMATE EXPONENTIAL UNDER- AND OVERFLOW LIMIT
C-----------------------------------------------------------------------
      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
      EXR = DCOS(ZR)
      EXI = DSIN(ZR)
      EY = 0.0D0
      TAY = DABS(ZI+ZI)
      IF (TAY.LT.ELIM) EY = DEXP(-TAY)
      IF (ZI.LT.0.0D0) GO TO 90
      C1R = EXR*EY
      C1I = EXI*EY
      C2R = EXR
      C2I = -EXI
   70 CONTINUE
      NZ = 0
      RTOL = 1.0D0/TOL
      ASCLE = D1MACH(1)*RTOL*1.0D+3
      DO 80 I=1,N
C       CY(I) = HCI*(C2*CWRK(I)-C1*CY(I))
C       STR = C1R*CYR(I) - C1I*CYI(I)
C       STI = C1R*CYI(I) + C1I*CYR(I)
C       STR = -STR + C2R*CWRKR(I) - C2I*CWRKI(I)
C       STI = -STI + C2R*CWRKI(I) + C2I*CWRKR(I)
C       CYR(I) = -STI*HCII
C       CYI(I) = STR*HCII
        AA = CWRKR(I)
        BB = CWRKI(I)
        ATOL = 1.0D0
        IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 75
          AA = AA*RTOL
          BB = BB*RTOL
          ATOL = TOL
   75   CONTINUE
        STR = (AA*C2R - BB*C2I)*ATOL
        STI = (AA*C2I + BB*C2R)*ATOL
        AA = CYR(I)
        BB = CYI(I)
        ATOL = 1.0D0
        IF (DMAX1(DABS(AA),DABS(BB)).GT.ASCLE) GO TO 85
          AA = AA*RTOL
          BB = BB*RTOL
          ATOL = TOL
   85   CONTINUE
        STR = STR - (AA*C1R - BB*C1I)*ATOL
        STI = STI - (AA*C1I + BB*C1R)*ATOL
        CYR(I) = -STI*HCII
        CYI(I) =  STR*HCII
        IF (STR.EQ.0.0D0 .AND. STI.EQ.0.0D0 .AND. EY.EQ.0.0D0) NZ = NZ
     *   + 1
   80 CONTINUE
      RETURN
   90 CONTINUE
      C1R = EXR
      C1I = EXI
      C2R = EXR*EY
      C2I = -EXI*EY
      GO TO 70
  170 CONTINUE
      NZ = 0
      RETURN
      END
      subroutine ZQCAI

c*********************************************************************72
c
C***DATE WRITTEN   830501   (YYMMDD)
C***REVISION DATE  890801, 930101   (YYMMDD)
C
C                *** A DOUBLE PRECISION ROUTINE ***
C
C     ZQCAI IS A QUICK CHECK ROUTINE FOR THE COMPLEX AIRY FUNCTIONS
C     GENERATED BY SUBROUTINES ZAIRY AND ZBIRY.
C
C     ZQCAI GENERATES AIRY FUNCTIONS AND THEIR DERIVATIVES FROM ZAIRY
C     AND ZBIRY AND CHECKS THEM AGAINST THE WRONSKIAN EVALUATION IN THE
C     REGION -PI/3 .LE. ARG(Z) .LE. PI/3:
C
C                 AI(Z)*BI'(Z)-AI'(Z)*BI(Z)=1/PI.
C
C     IN THE REMAINDER OF THE CUT PLANE, THE IDENTITIES
C
C              AI(Z)  = SQRT(-Z)*( J(-1/3,ZR) + J(1/3,ZR) )/3
C
C              AI'(Z) =        Z*( J(-2/3,ZR) - J(2/3,ZR) )/3
C
C       BI(Z)  =   I*SQRT(-Z/3)*( C1*H(1/3,1,ZR) - C2*H(1/3,2,ZR) )/2
C
C       BI'(Z) = I*(-Z)/SQRT(3)*( C2*H(2/3,1,ZR) - C1*H(2/3,2,ZR) )/2
C
C     ARE CHECKED WHERE ZR = (2/3)(-Z)**(3/2) WITH C1 = EXP(PI*I/6),
C     C2 = CONJG(C1) AND I**2 = -1.
C
C     THE PARAMETER MQC CAN HAVE VALUES 1 (THE DEFAULT) FOR A FASTER,
C     LESS DEFINITIVE TEST OR 2 FOR A SLOWER, MORE DEFINITIVE TEST.
C
C     MACHINE CONSTANTS ARE DEFINED IN FUNCTIONS I1MACH, R1MACH, AND
C     D1MACH. THESE MUST BE SELECTED BY THE USER OR SET ACCORDING TO
C     PROLOGUE INSTRUCTIONS.
C
C     COMPLEX CA, CAV, CHI, CI, CONA, CONB, CONC, COND, CON1, CON2, CV,
C    *  CW, CY, W, Y, YY, Z, ZR, ZW, SC
      EXTERNAL ZABS
      DOUBLE PRECISION AA, ALIM, ATOL, CAR, CAI, CVR, CVI,
     * CON1I, CON1R, CON2I, CON2R, CONAI, CONAR, CONBI, CONBR, CONCI,
     * CONCR, CONDI, CONDR, C13, C23, C43, CAVR, CAVI, CIR, CII, CHIR,
     * CHII, CWI, CWR, CYI, CYR, CT, DIG, ELIM, EPS, ER, ERTOL, FNUL,
     * FPI, HPI, PI, PTR, R, RL, RPI, R1M4, R1M5, SPI, STI, STR, ST,
     * T, TOL, RM, WI, WR, YI, YR, YYR, YYI, ZI, ZR, ZRI, ZRR, ZWR, ZWI,
     * D1MACH, ZABS, PI3, SLAK, FILM, AB, TS, SCR, SCI
      INTEGER I, ICASE, IL, IR, IRB, IRSET, IT, ITL, K, KEPS, KODE, J,
     * JB, JL, K1, K2, LFLG, LUN, NZ, IERR, I1MACH, MQC, ICL, KDO
      DIMENSION ER(5), T(20), YR(20), YI(20), YYR(20), YYI(20), WR(20),
     * WI(20), KEPS(20), KDO(20)

      PARAMETER (MQC=1)

C-----------------------------------------------------------------------
C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
C-----------------------------------------------------------------------
      R1M4 = D1MACH(4)
      TOL = MAX(R1M4,1.0D-18)
      AA = -DLOG10(R1M4)
      K1 = I1MACH(15)
      K2 = I1MACH(16)
      R1M5 = D1MACH(5)
      K = MIN(ABS(K1),ABS(K2))
      ELIM = 2.303D0*(DBLE(FLOAT(K))*R1M5-3.0D0)
      AB = AA*2.303D0
      ALIM = ELIM + MAX(-AB,-41.45D0)
      DIG = MIN(AA,18.0D0)
      FNUL = 10.0D0 + 6.0D0*(DIG-3.0D0)
      RL = 1.2D0*DIG + 3.0D0
      SLAK = 3.0D0+4.0D0*(-DLOG10(TOL)-7.0D0)/11.0D0
      SLAK = MAX(SLAK,3.0D0)
      ERTOL = TOL*10.0D0**SLAK
      RM = 0.5D0*(ALIM + ELIM)
      RM=MIN(RM,200.0D0)
      RM=MAX(RM,RL+10.0D0)
C-----------------------------------------------------------------------
      write ( *,99999)
99999 FORMAT (' QUICK CHECK ROUTINE FOR THE AIRY FUNCTIONS FROM ZAIRY AN
     *D ZBIRY'/)
      write ( *,99998)
99998 FORMAT (' PARAMETERS TOL,ELIM,ALIM,RL,FNUL,DIG')
      write ( *,99997) TOL, ELIM, ALIM, RL, FNUL, DIG
99997 FORMAT (6D12.4/)
      ATOL = 100.0D0*TOL
      FPI = DATAN(1.0D0)
      HPI = FPI + FPI
      PI = HPI + HPI
      RPI = 1.0D0/PI
      SPI = PI/6.0D0
      CON1R = COS(SPI)
      CON1I = SIN(SPI)
      CON2R = CON1R
      CON2I = -CON1I
      PI3 = SPI+SPI
      C13 = 1.0D0/3.0D0
      C23 = C13+C13
      C43 = C23+C23
      CAVR = SQRT(C13)
      CAVI = 0.0D0
      CHIR = 0.0D0
      CHII = 0.5D0
      CIR = 0.0D0
      CII = 1.0D0
      write ( *,99996) MQC
99996 FORMAT (/' CHECKS IN THE (Z,FNU) SPACE WITH MQC = ',I2/)
C-----------------------------------------------------------------------
C     TEST VALUES OF Z IN -PI.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     KDO(K), K=1,IL  DETERMINES WHICH OF THE IL ANGLES IN -PI TO PI
C     ARE USE TO COMPUTE VALUES OF Z
C       KDO(K) = 0  MEANS THAT THE INDEX K WILL BE USED FOR ONE OR TWO
C                   VALUES OF Z, DEPENDING ON THE CHOICE OF KEPS(K)
C              = 1  MEANS THAT THE INDEX K AND THE CORRESPONDING ANGLE
C                   WILL BE SKIPPED
C     KEPS(K), K=1,IL DETERMINES WHICH OF THE ANGLES GET INCREMENTED
C     UP AND DOWN TO PUT VALUES OF Z IN REGIONS WHERE DIFFERENT
C     FORMULAE ARE USED.
C       KEPS(K) =0  MEANS THAT THE ANGLE WILL BE USED WITHOUT CHANGE
C               =1  MEANS THAT THE ANGLE WILL BE INCREMENTED UP AND
C                   DOWN BY EPS
C     THE ANGLES TO BE USED ARE STORED IN THE T(I) ARRAY, I=1,ITL
C-----------------------------------------------------------------------
      IF (MQC.NE.2) THEN
        ICL=1
        IL=5
        DO 5 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    5   CONTINUE
      ELSE
        ICL=2
        IL=7
        DO 6 I=1,IL
          KDO(I)=0
          KEPS(I)=0
    6   CONTINUE
        KEPS(2)=1
        KEPS(3)=1
        KEPS(5)=1
        KEPS(6)=1
      ENDIF
      I = 2
      EPS=0.01D0
      FILM=DBLE(FLOAT(IL-1))
      T(1) = -PI + EPS
      DO 30 K=2,IL
        IF(KDO(K).EQ.0) THEN
          T(I) = PI*DBLE(FLOAT(-IL+2*K-1))/FILM
          IF (KEPS(K).EQ.0) GO TO 20
          TS=T(I)
          T(I) = TS - EPS
          I = I + 1
          T(I) = TS + EPS
   20     CONTINUE
          I = I + 1
        ENDIF
   30 CONTINUE
      ITL = I - 1
      LFLG = 0
      DO 180 ICASE=1,ICL
        DO 170 KODE=1,2
          DO 160 IRSET=1,3
            IRB = MIN(IRSET,2)
            DO 150 IR=IRB,4
              GO TO (40, 50, 60), IRSET
   40         CONTINUE
              R =(0.2D0*DBLE(FLOAT(4-IR))+2.0D0*DBLE(FLOAT(IR-1)))/3.0D0
              GO TO 70
   50         CONTINUE
              R = (2.0D0*DBLE(FLOAT(4-IR))+RL*DBLE(FLOAT(IR-1)))/3.0D0
              GO TO 70
   60         CONTINUE
              R = (RL*DBLE(FLOAT(4-IR))+RM*DBLE(FLOAT(IR-1)))/3.0D0
   70         CONTINUE
              DO 140 IT=1,ITL
                CT = COS(T(IT))
                ST = SIN(T(IT))
                IF (DABS(CT).LT.ATOL) CT = 0.0D0
                IF (DABS(ST).LT.ATOL) ST = 0.0D0
                ZR = R*CT
                ZI = R*ST
                IF(ABS(T(IT)).LE.PI3) THEN
C-----------------------------------------------------------------------
C     WRONSKIAN CHECK IN -PI/3.LT.ARG(Z).LT.PI/3, TEST #1
C-----------------------------------------------------------------------
                  CALL ZAIRY(ZR, ZI, 0, KODE, YR(1), YI(1), NZ, IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL ZAIRY(ZR, ZI, 1, KODE, YR(2), YI(2), NZ, IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL ZBIRY(ZR, ZI, 0, KODE, WR(1), WI(1), IERR)
                  CALL ZBIRY(ZR, ZI, 1, KODE, WR(2), WI(2), IERR)
                  CWR = YR(1)*WR(2)-YI(1)*WI(2)
                  CWI = YR(1)*WI(2)+YI(1)*WR(2)
                  CYR = YR(2)*WR(1)-YI(2)*WI(1)
                  CYI = YR(2)*WI(1)+YI(2)*WR(1)
                  CVR = RPI
                  CVI = 0.0D0
                  IF (KODE.EQ.2) THEN
                    CALL ZSQRT(ZR,ZI,CAR,CAI)
                    ZRR = (ZR*CAR-ZI*CAI)*C23
                    ZRI = (ZR*CAI+ZI*CAR)*C23
                    AA=ABS(ZRR)
                    CAR = ZRR - AA
                    CAI = ZRI
                    CALL ZEXP(CAR,CAI,STR,STI)
                    PTR = STR*CVR-STR*CVI
                    CVI = STR*CVI+STI*CVR
                    CVR = PTR
                  ENDIF
                  CYR = CWR - CYR - CVR
                  CYI = CWI - CYI - CVI
                  ER(1) = ZABS(CYR,CYI)/ZABS(CVR,CVI)
                  JB = 1
                  JL = 1
                ELSE
C-----------------------------------------------------------------------
C     CHECKS IN -PI.LT.ARG(Z).LT.-PI/3 AND PI/3.LT.ARG(Z).LE.PI
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     CHECK AI    TEST #2
C-----------------------------------------------------------------------
                  CALL ZAIRY(ZR, ZI, 0, KODE, YR(2), YI(2), NZ, IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  ZRR = -ZR
                  ZRI = -ZI
                  CALL ZSQRT(ZRR,ZRI,CVR,CVI)
                  PTR = (ZRR*CVR-ZRI*CVI)*C23
                  ZRI = (ZRR*CVI+ZRI*CVR)*C23
                  ZRR = PTR
                  CALL ZBESJ(ZRR,ZRI,C23,KODE,2,YYR,YYI,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  STR = YYR(1)*C43
                  STI = YYI(1)*C43
                  CALL ZDIV(STR,STI,ZRR,ZRI,STR,STI)
                  CYR = STR - YYR(2)
                  CYI = STI - YYI(2)
                  CAR = YYR(1)
                  CAI = YYI(1)
                  CALL ZBESJ(ZRR,ZRI,C13,KODE,2,YYR,YYI,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  IF (KODE.EQ.2) THEN
                    AB = ABS(ZRI)
                    CALL ZSQRT(ZR,ZI,CWR,CWI)
                    ZWR = (ZR*CWR-ZI*CWI)*C23
                    ZWI = (ZR*CWI+ZI*CWR)*C23
                    CWR = ZWR+AB
                    CWI = ZWI
                    CALL ZEXP(CWR,CWI,STR,STI)
                    CWR = STR
                    CWI = STI
                    STR = YYR(1)*CWR-YYI(1)*CWI
                    YYI(1) = YYR(1)*CWI+YYI(1)*CWR
                    YYR(1) = STR
                    STR = YYR(2)*CWR-YYI(2)*CWI
                    YYI(2) = YYR(2)*CWI+YYI(2)*CWR
                    YYR(2) = STR
                    STR = CYR*CWR-CYI*CWI
                    CYI = CYR*CWI+CYI*CWR
                    CYR = STR
                    STR = CAR*CWR-CAI*CWI
                    CAI = CAR*CWI+CAI*CWR
                    CAR = STR
                    SCR = CWR
                    SCI = CWI
                  ENDIF
                  CWR = CVR*C13
                  CWI = CVI*C13
                  WR(2) = CWR*(YYR(1)+CYR)-CWI*(YYI(1)+CYI)
                  WI(2) = CWR*(YYI(1)+CYI)+CWI*(YYR(1)+CYR)
                  STR = YR(2)-WR(2)
                  STI = YI(2)-WI(2)
                  ER(2) = ZABS(STR,STI)
                  IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN
                    ER(2) = ER(2)/ZABS(YR(2),YI(2))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(2) = ER(2)/ZABS(SCR,SCI)
                    ENDIF
                  ENDIF
C-----------------------------------------------------------------------
C     CHECK AI'   TEST #3
C-----------------------------------------------------------------------
                  STR = YYR(1)*C23
                  STI = YYI(1)*C23
                  CALL ZDIV(STR,STI,ZRR,ZRI,CYR,CYI)
                  CYR = CYR-YYR(2)-CAR
                  CYI = CYI-YYI(2)-CAI
                  WR(3) = (ZR*CYR-ZI*CYI)*C13
                  WI(3) = (ZR*CYI+ZI*CYR)*C13
                  CALL ZAIRY(ZR, ZI, 1, KODE, YR(3), YI(3), NZ, IERR)
                  STR = YR(3)-WR(3)
                  STI = YI(3)-WI(3)
                  ER(3) = ZABS(STR,STI)
                  IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN
                    ER(3) = ER(3)/ZABS(YR(3),YI(3))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(3) = ER(3)/ZABS(SCR,SCI)
                    ENDIF
                  ENDIF
C-----------------------------------------------------------------------
C     CHECK BI    TEST #4
C-----------------------------------------------------------------------
                  CALL ZBESH(ZRR,ZRI,C13,KODE,1,1,YR,YI,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL ZBESH(ZRR,ZRI,C13,KODE,2,1,YYR,YYI,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CONAR = CON1R
                  CONAI = CON1I
                  CONBR = CON2R
                  CONBI = CON2I
                  CONCR = CON2R
                  CONCI = CON2I
                  CONDR = CON1R
                  CONDI = CON1I
                  IF (KODE.EQ.2) THEN
                    AA = ABS(ZWR)
                    ZWR = CIR*ZRR-CII*ZRI-AA
                    ZWI = CIR*ZRI+CII*ZRR
                    CALL ZEXP(ZWR,ZWI,CWR,CWI)
                    STR = CONAR*CWR-CONAI*CWI
                    CONAI =CONAR*CWI+CONAI*CWR
                    CONAR = STR
                    STR = CONCR*CWR-CONCI*CWI
                    CONCI =CONCR*CWI+CONCI*CWR
                    CONCR = STR
                    ZWR = -(CIR*ZRR-CII*ZRI)-AA
                    ZWI = -(CIR*ZRI+CII*ZRR)
                    CALL ZEXP(ZWR,ZWI,CWR,CWI)
                    STR = CONBR*CWR-CONBI*CWI
                    CONBI =CONBR*CWI+CONBI*CWR
                    CONBR = STR
                    STR = CONDR*CWR-CONDI*CWI
                    CONDI =CONDR*CWI+CONDI*CWR
                    CONDR = STR
                    SCR = CWR
                    SCI = CWI
                  ENDIF
                  CWR = CONAR*YR(1)-CONAI*YI(1)
                  CWI = CONAR*YI(1)+CONAI*YR(1)
                  CWR = CWR - (CONBR*YYR(1)-CONBI*YYI(1))
                  CWI = CWI - (CONBR*YYI(1)+CONBI*YYR(1))
                  STR = CVR*CAVR-CVI*CAVI
                  STI = CVR*CAVI+CVI*CAVR
                  PTR = STR*CWR-STI*CWI
                  CWI = STR*CWI+STI*CWR
                  CWR = PTR
                  WR(4) = CWR*CHIR-CWI*CHII
                  WI(4) = CWR*CHII+CWI*CHIR
                  CALL ZBIRY(ZR,ZI,0,KODE,YR(4),YI(4),IERR)
                  STR = YR(4)-WR(4)
                  STI = YI(4)-WI(4)
                  ER(4) = ZABS(STR,STI)
                  IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN
                    ER(4) = ER(4)/ZABS(YR(4),YI(4))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(4) = ER(4)/ZABS(SCR,SCI)
                    ENDIF
                  ENDIF
C-----------------------------------------------------------------------
C     CHECK BI'   TEST #5
C-----------------------------------------------------------------------
                  CALL ZBESH(ZRR,ZRI,C23,KODE,1,1,YR,YI,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CALL ZBESH(ZRR,ZRI,C23,KODE,2,1,YYR,YYI,NZ,IERR)
                  IF (NZ.NE.0.OR.IERR.NE.0) GO TO 140
                  CWR = CONCR*YR(1)-CONCI*YI(1)
                  CWI = CONCR*YI(1)+CONCI*YR(1)
                  CWR = CWR - (CONDR*YYR(1)-CONDI*YYI(1))
                  CWI = CWI - (CONDR*YYI(1)+CONDI*YYR(1))
                  STR = -(ZR*CAVR-ZI*CAVI)
                  STI = -(ZR*CAVI+ZI*CAVR)
                  PTR = STR*CWR-STI*CWI
                  CWI = STR*CWI+STI*CWR
                  CWR = PTR
                  WR(5) = CWR*CHIR-CWI*CHII
                  WI(5) = CWR*CHII+CWI*CHIR
                  CALL ZBIRY(ZR,ZI,1,KODE,YR(5),YI(5),IERR)
                  STR = YR(5)-WR(5)
                  STI = YI(5)-WI(5)
                  ER(5) = ZABS(STR,STI)
                  IF (ZI.NE.0.0D0.OR.ZR.GE.0.0D0) THEN
                    ER(5) = ER(5)/ZABS(YR(5),YI(5))
                  ELSE
                    IF (KODE.EQ.2) THEN
                      ER(5) = ER(5)/ZABS(SCR,SCI)
                    ENDIF
                  ENDIF
                  JB = 2
                  JL = 5
                ENDIF
                DO 190 J=JB,JL
                IF (ER(J).LT.ERTOL) GO TO 190
                IF (LFLG.EQ.1) GO TO 130
                write ( *,99995) ERTOL
99995           FORMAT (/' CASES WHICH VIOLATE THE RELATIVE ERROR TEST W
     *ITH ERTOL =', E12.4/)
                write ( *,99994)
99994           FORMAT (/' OUTPUT FORMAT'/' KODE,IR,IT,IRSET,ICASE')
                write ( *,99993)
99993           FORMAT (' ER'/' I, Z, Y(I), W(I), ON THE I-TH TEST, I=1,
     *5'/)
                LFLG = 1
  130           CONTINUE
                write ( *,99992) KODE, IR, IT, IRSET, ICASE
99992           FORMAT (5I5)
                write ( *,99991) ER(J)
                write ( *,99989) J, ZR, ZI, YR(J), YI(J), WR(J), WI(J)
99991           FORMAT (D12.4)
99989           FORMAT (I5,6D12.4)
  190           CONTINUE
  140         CONTINUE
  150       CONTINUE
  160     CONTINUE
  170   CONTINUE
  180 CONTINUE
      IF (LFLG.EQ.0) write ( *,99990)
99990 FORMAT (/' QUICK CHECKS OK'/)
      return
      END
